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SECTION  1 
INTRODUCTION 

Satellite  communications  systems  that  utilize  transionospheric  propagation  links 
may  be  subject  to  severe  performance  degradation  when  the  ionosphere  is  highly 
disturbed  by  high  altitude  nuclear  explosions  1  Arendt  and  Soichcr  1964;  King  and 
Fleming  1980]  or  by  chemical  releases  \ Davis  et  al.  1974;  Wolcott  et  al.  1978).  During 
these  events,  the  increased  electron  concentrations  and  the  irregular  structure  of  the 
ionization  can  lead  to  intense  Rayleigh  signal  scintillation  at  the  radio  frequencies  (RF) 
used  for  satellite  communication  links  and  space  radars. 

Under  severe  scintillation  conditions,  the  signal  incident  at  the  receiver  can  vary 
randomly  in  amplitude,  phase,  time-of-arrival,  and  angle-of-arrival.  If  all  frequency 
components  of  the  signal  vary  essentially  identically  w'ith  time,  the  propagation  channel 
is  referred  to  as  nonselective  or  flat  fading.  When  the  scintillations  exhibit  statistical 
decorrelation  at  different  frequencies  within  the  signal  bandwidth,  the  channel  is 
referred  to  as  frequency  selective.  Frequency  selective  scintillations  are  therefore 
encountered  when  the  signal  bandwidth  exceeds  the  frequency  selective  bandwidth  of 
the  channel.  When  the  scintillations  exhibit  statistical  decorrelation  across  the  face  of 
an  aperture  antenna,  the  channel  may  also  be  referred  to  as  spatially  selective. 
Spatially  selective  scintillations  are  therefore  encountered  when  the  antenna  aperture 
size  exceeds  the  decorrelation  distance  of  the  incident  signal. 

Under  conditions  where  the  signal  is  spatially  selective,  the  antenna  beamwidth  is 
smaller  than  the  angle-of-arrival  fluctuations  and  the  effect  of  the  antenna  is  to 
attenuate  the  incident  signal  that  is  arriving  at  off-boresight  angles.  In  the  spatial 
domain,  the  incident  electric  field  is  somewhat  decorrelated  across  the  face  of  the 
antenna.  The  induced  voltages  in  the  antenna  then  do  not  add  coherently  as  they  would 
for  an  incident  plane  wave,  resulting  in  a  loss  in  the  gain  of  the  antenna.  Because  of 
this  angular  filtering  or  spatial  selectivity,  the  second-order  statistics  of  the  signal  at  the 
output  of  the  antenna  will  be  different  than  those  of  the  incident  signal. 

The  effects  of  antennas  on  signals  that  have  propagated  through  randomly 
ionized  media  have  been  reported  by  Wittwer  [1982],  Knepp  [1985],  and  Dana  ]  1986]. 
This  report  is  an  extension  of  the  latter  reference,  and  its  purpose  is  to  describe  the 
general  model  for  temporal  fluctuations  that  has  recently  been  developed  by  Dr.  Leon 
A.  Wittwer  of  the  Defense  Nuclear  Agency.  A  review  of  the  basic  theory  of  RF 
propagation  through  random  media  is  presented,  and  the  channel  simulation  technique 
for  the  general  model  is  described. 

The  starting  point  is  the  generalized  power  spectral  density  (GPSD).  The  first 
part  of  this  report  is  a  review  of  the  derivation  of  the  GPSD.  The  intent  of  this  review 
is  to  give  the  reader  an  understanding  of  the  underlying  physics  that  are  contained  in 
the  GPSD  and  an  understanding  of  the  assumptions  used  to  calculate  die  GPSD.  The 
first  part  of  this  review  follows  Tatarskii  1 1 97 1 ,  §64-65].  The  discussion  of  the 
general  model  is  new  in  this  report. 

The  derivation  of  the  GPSD  starts  with  Maxwell’s  equations  from  which  the 
parabolic  wave  equation  is  derived.  The  parabolic  wave  equation  can  be  solved  to  give 
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the  received  electric  field  for  a  specific  electron  density  distribution  in  the  ionosphere. 
However,  the  electron  density  distribution  is  a  random  process  so  the  received  electric 
field  is  also  a  random  process.  The  parabolic  wave  equation  is  therefore  used  to  derive 
an  equation  for  the  two-position,  tw'o-frequency,  two-time  mutual  coherence  function 
of  the  electric  field,  r(8r,8co,8t).  The  solution  of  the  differential  equation  for  T, 
which  is  also  a  solution  of  Maxwell's  equations,  then  provides  a  description  of  the 
second-order  statistics  of  the  received  electric  field.  The  Fourier  transform  of  the 
mutual  coherence  funcfion  is  the  GPSD  of  the  received  signal. 

Once  the  GPSD  is  obtained  for  the  general  model,  it  is  used  to  compute  the  mean 
power,  decorrelation  time,  and  frequency  selective  bandwidth  of  the  signal  out  of 
anisotropic  antennas  with  arbitrary  pointing  angles  relative  to  the  line-of-sight.  These 
results  are  taken  from  Frasier  [1988].  Several  examples  are  given  that  illustrate  the 
general  model  and  the  effects  of  antenna  pointing. 

In  Section  4  an  analytical/numerical  technique  is  described  that  is  used  to 
generate  realizations  of  the  impulse  response  function  of  the  signal  after  propagation 
through  randomly  ionized  media  and  reception  by  multiple  antennas.  The  statistical 
realizations  of  the  signal  at  the  outputs  of  multiple  antennas  are  assumed  to  have 
Rayleigh  amplitude  statistics  and  are  therefore  valid  under  strong-scattering  conditions. 
The  spatial  and  frequency  correlation  properties  of  the  realizations  are  given  by  the 
mutual  coherence  function.  These  realizations  of  the  impulse  response  function  are 
then  used  to  construct  the  received  signal  and  may  be  used  to  exercise  simulations  of 
transionospheric  communications  links  or  space  radars. 


SECTION  2 
THEORY 


The  starting  point  for  this  discussion  of  the  general  model  is  the  generalized 
power  spectral  density  (GPSD).  This  section  presents  a  review  of  the  derivation  of  the 
GPSD  and  discusses  the  physics  that  are  contained  in  this  important  function. 

In  deriving  the  GPSD,  two  key  approximations  are  usually  made  about  the 
spatial  and  temporal  electron  density  fluctuations  that  cause  the  scattering  in  the 
ionosphere.  The  first  of  these  is  the  delta-layer  approximation  which  requires  that  the 
scattering  occur  in  an  infinitesimally  thin  layer  normal  to  the  line-of-sight.  This 
approximation  has  been  relaxed  in  the  calculations  of  Witiwer  [1982]  and  Knepp  [1983] 
and  has  been  found  to  result  in  small  errors  in  the  GPSD  provided  that  the  propagation 
parameters  (frequency  selective  bandwidth,  decorrelation  time  and  decorrelation 
distance)  that  characterize  the  channel  are  properly  specified.  The  delta-layer 
approximation  is  not,  in  general,  adequate  to  calculate  the  propagation  parameters. 

The  second  approximation  that  is  usually  made  is  Taylor's  frozen-in  hypothesis 
which  treats  the  ionization  fluctuations  or  striations  as  rigid  “frozen-in”  structures  that 
drift  across  the  line-of-sight.  Under  this  model  there  is  strong  coupling  between  the 
spatial  and  temporal  variations  of  the  random  electric  field  that  is  incident  at  the  plane 
of  the  receiver.  Wittwer  [1988]  has  recently  proposed  a  “general  model”  that  smoothly 
varies  between  the  frozen-in  model  and  a  turbulent  model  where  the  temporal  and 
spatial  ionization  fluctuations  are  uncorrelated.  The  GPSD  for  the  general  model  will 
be  calculated  in  this  section. 

An  analytic  solution  is  obtained  in  this  section  for  the  two-position,  two- 
frequency,  two-time  mutual  coherence  function  r(5r,8co,5t)  of  the  complex  electric 
field  incident  on  the  plane  of  the  receiving  antenna.  This  solution  is  valid  for  arbitrary 
line-of-sight  geometries  relative  to  the  ionization  structures  in  the  ionosphere  that  cause 
the  scattering  of  the  RF  wave.  The  mutual  coherence  function  then  provides  the  basis 
for  the  antenna  aperture  effects  calculations  and  for  the  statistical  signal  generation 
techniques  discussed  in  subsequent  sections  of  this  report.  The  Fourier  transform  of 
the  mutual  coherence  function  is  the  GPSD  of  the  received  signal. 

2.1  PARABOLIC  WAVE  EQUATION. 

Consider  a  monochromatic  spherical  wave  with  an  electric  field  £(r,co,t)  which 
is  a  function  of  position  r,  carrier  radian  frequency  co,  and  time  t.  The  wave 
originates  from  a  transmitter  located  at  r  =  (0,0, -z^)  and  propagates  in  free  space  in  the 
positive  z  direction  until  it  is  incident  on  an  irregularly  ionized  layer  which  extends 
from  0  <  z  <  L  and  is  infinite  in  the  x-y  plane.  After  emerging  from  the  layer  at 
z  =  L,  the  wave  propagates  in  free  space  to  a  receiver  located  at  r  =  (0,0, zr).  This 
geometry  is  shown  in  Figure  1 . 
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Figure  1.  Propagation  geometry. 


The  propagation  of  the  wave  is  governed  by  Maxwell's  equations: 

***  +  ''%-  o  VE  =  0 

VxH'c¥  =  0  VH  =  0 


(1) 


where  c  is  the  speed  of  light  in  vacuum,  H  is  the  magnetic  field,  and  e  is  the  dielectric 
constant. 

The  dielectric  constant  of  the  propagation  medium  undergoes  random 
fluctuations  with  a  characteristic  frequency  which  is  assumed  to  be  small  when 
compared  to  the  carrier  frequency  of  the  wave.  With  this  assumption,  the  electric  and 
magnetic  fields  may  be  written  as  the  product  of  slowly-varying  complex  envelopes, 
denoted  E  and  H,  times  exp(io)t): 

£(r,co,t)  =  E(r,(o,t)  eiwt  (2) 

H(r,o),t)  =  H(r,co,t)  eIOM  . 

Inserting  these  into  Maxwell's  equations  gives 

V  x  E  +  ikH  =0  (3a) 

V  x  H  -  iekE  =  0  (3b) 
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where  k  =  co/c  is  the  wave  number  of  the  carrier.  After  applying  the  curl  operator  to 
Equation  3a  and  substituting  Equation  3b  for  the  V  x  H  term,  the  equation  for  E 
becomes 

V  x  V  x  E  -  e  k;  E  =  0  .  (4) 

The  vector  identity  V  x  V  x  E  =  V(V-E)  -  V;E  reduces  the  curl  curl  term  in  Equation 
4  with  the  result 

V  (V  E)  -  V2E  -  e  k2  E  =  0  .  (5) 

The  V  E  term  is  reduced  by  expanding  the  divergence  equation  for  E: 


V  cE  =  e  V  E  +  E  V  e  =  0 


(6) 


or 


V  E  =  -  E  -V  (In  e)  . 


(7) 


The  wave  equation  for  the  complex  envelope  of  the  electric  field  then  becomes 
V:E  +  e  k2  E  +  V  [E  V  (In  e)]  =  0  . 


(8) 


The  dielectric  constant  e  in  a  plasma  at  radio  frequencies  is  given  approximately 
by 


*» 


(9) 


where  the  plasma  frequency  is 

co n  =  4rc  re  c2  ne(r,t)  .  (10) 

The  quantity  re  is  the  classical  electron  radius  (re  =  2.8179x10  15  m)  and  ne(r,t)  is  the 
free  electron  density  as  a  function  of  position  and  time.  Equation  9  is  valid  when  the 
carrier  frequency  is  large  compared  to  the  plasma  frequency.  The  free  electron 
density  is  a  random  variable  that  will  be  represented  as  a  mean  value  plus  a  random 
variation: 


ne(r,t)  =  (ne>  +  Ane(r,t)  .  (11) 

The  electron  density  fluctuation  Ane(r,t)  is  assumed  to  be  a  zero  mean  random  process. 
The  term  e  k2  in  the  wave  equation  may  now  be  rewritten  as 

tk‘  =  (k)2  (1  -  e, )  (12) 


where 
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4n:rcc:(ne)  ] 
or  J 


and 


"Anc(r,tf 

<«p> 

U  <ne>  J 

_co2  -  (orp)j 

(13) 


(14) 


The  quantity  <cn^}  is  the  plasma  frequency  evaluated  at  the  mean  electron  density. 

The  magnitude  of  the  gradient  term  V(In  e)  in  the  wave  equation  may  be 
estimated  as  follows: 


lV(ln  e)l 


1  - 


(Oj, 

to2 


47rre(nc}c2 

co2 


<<0p)  1 

CO2  U 

(15) 


where  L0  is  the  scale  size  of  the  electron  density  fluctuations.  As  long  as  L0  »  X, 
where  X  is  the  RF  wavelength,  the  term  V[EV(ln  £)]  is  small  compared  to  ek"E  and 
may  be  ignored.  The  steps  that  follow  are  therefore  only  valid  when  the  scale  size  of 
the  electron  density  fluctuations  is  large  compared  to  the  carrier  wavelength.  With  this 
restriction,  the  wave  equation  takes  the  Helmholtz  form: 

V2E  +  (k)2  (1  -  e,)E  =  0  .  (16) 


Now  consider  the  complex  components  of  the  electric  field  and  let 

E(r,co,t)  =  U(r,co,t)  exp  [- J  <k(z')>  dz' ]  .  (17) 


This  scalar  equation  for  E  may  be  used  because  it  is  usual  for  trans-ionospheric  RF 
links  to  be  circularly  polarized.  It  is  therefore  not  necessary  to  carry  out  separate 
calculations  for  each  polarization  state.  The  exponential  term  in  Equation  17 
represents  the  dispersive  effects  of  the  smooth  plasma  and  will  be  discussed  later  in  this 
section.  The  electric  field  envelope  U  contains  the  diffractive  effects  that  are  of 
interest  under  strong-scattering  conditions.  Substituting  Equation  17  for  E(r,co,t)  in 
the  wave  equation  gives  the  following  differential  equation  for  U: 

VJu+0-2i(k>^-<k)2E,U  =  0  (18) 


where 


2  _  ai  ai 

1  “  5x2  +  dy2 


(19) 


The  complex  amplitude  U  varies  as  the  electron  density  fluctuations.  The  second 
derivative  d2U/dz2  is  then  the  order  of  U/Lq.  On  the  other  hand,  the  term  (k)  dU/dz 
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varies  as  U/XL{).  As  long  as  X,  «  L0,  the  second  derivative  is  small  compared  to  the 
first  derivative  and  may  be  ignored.  This  is  equivalent  to  neglecting  reflected  rays  and 
is  called  the  “parabolic”  approximation.  The  parabolic  wave  equation  is  then 

VlU-2i(k)^-(k)2e,U  =  0  .  (20) 

It  will  be  seen  that  this  parabolic  or  small-angle  scattering  approximation  is  robust  in 
that  it  degrades  gracefully  as  the  scattering  angles  get  large.  The  source  term  £)U  in 
the  parabolic  wave  equation  is  a  function  of  frequency  and  the  electron  density 
fluctuations.  Different  frequencies  within  a  signal  bandwidth  may  therefore  propagate 
differently  through  the  same  ionization  structure.  When  this  happens,  the  propagation 
channel  is  said  to  be  frequency  selective. 

2.2  TRANSPORT  EQUATION. 

A  partial  differential  equation  for  the  two-position,  two-frequency,  two-time 
mutual  coherence  function  is  derived  from  the  parabolic  wave  equation  in  this  section. 
This  transport  equation  is  derived  using  the  Novikov  theorem  which  requires  that  the 
electron  density  fluctuations  be  normally  distributed.  However,  Lee  and  Jokipii 
[  1 975a]  give  an  alternative  derivation  that  relaxes  this  assumption. 

2.2.1  First  Form  of  the  Transport  Equation. 

The  two-position,  two- frequency,  two-time  mutual  coherence  function  is  defined 
in  a  plane  normal  to  the  line-of-sight  as 

T  =  <U(p,,z,©i,ti)  U*(p2,z,(02,t2)>  (21) 


where  p  is  a  two-dimensional  position  vector  in  the  plane  normal  to  the  line-of-sight. 

In  order  to  obtain  an  equation  for  T,  the  parabolic  wave  equation  for 
Ufp^z.coj.tj)  is  multiplied  by  U*(p2,z,a>2,t2).  This  results  in  the  following  equation: 


ViiU(p1,z,0)1,t1)U*(p2,z,0)2,t2)  -  2i 


5U(pltz,(0i,ti) 

dz 


U*(p2,z,co2,t2) 


-  kj  e^p^z.coj.tj)  U(p1,z,o)1,t1)  U*(p2,z,oo2,t2)  =  0  (22a) 

where  kj  is  given  by  Equation  13  evaluated  at  frequency  tOj,  (Pj,z,00j,tj)  is  given  by 
Equation  14,  and  the  Laplacian  is  given  by  Equation  19  evaluated  at  pj.  A  similar 
equation  can  be  written  down  by  interchanging  the  subscripts  1  and  2  and  by  taking  the 
complex  conjugate  with  the  result: 


1 

k2 


VJ%U*(p2,z,co2,t2)U(Pi,z,(Oi,t1)  +  2i 


dU*(p2,z,(02,t2) 

dz 


U  (p  i .  z,  co ,  ,t  i ) 


-  k2  £i  (p2>z>w2’l2)  U*(p2,Z,C02,t2)  U(Pi,Z,C0],ti )  -  0  .  (22b) 

Upon  subtracting  Equation  22b  from  Equation  22a  and  taking  the  expectation  value, 
the  equation  for  T  is 
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(23) 


_  i  n2  r  ..  or 
11  r  -  ^  v±2  r  -  2i  — 


k2  T  12  1  '  dz 
(CllPwZ.Wj.tj)  U(pi,Z,Ci)1,t1)  U*(p2,z,w2,t2)> 


+  k2  <Ei(p2,z,0)2,t2)  U*(p2,z,co2,t2)  U(p1,z,co1,t1)>  =  0  . 

The  expectation  of  the  two  source  terms  in  this  equation  must  be  carefully 
evaluated.  They  involve  the  product  of  UU*  and  £]  where  is  proportional  to  the 
fluctuations  in  the  electron  density.  However,  the  electric  field  complex  envelope  U  is 
also  a  function  of  the  electron  density  fluctuations  that  are  encountered  along  the 
propagation  path.  Therefore  U  and  e,  are  correlated. 

2.2.2  Novikov  Theorem. 

The  Novikov  theorem  is  used  to  evaluate  the  source  terms  in  Equation  23.  This 
theorem  is  proven  in  Tatarskii  [1971,  §65 [  and  Ishimaru  [1978,  pp.  457-458).  The 
theorem  states  that 


<f,(R)f2(fi)> 


<f  t  (R  )f ,  (R  ')> 


dnR’ 


(24) 


where  fj(R)  is  a  zero-mean,  normally-distributed  random  function  of  the  n- 
dimensional  vector  R,  f2( f j )  is  a  function  of  f,(  and  5f2/8fj  is  a  functional  derivative. 
In  applying  this  theorem,  f,  =  £,  and  f2  =  UU*.  The  theorem  is  proven  by  expanding 
f2(fj)  in  a  Taylor  series. 

2.2.3  Source  Terms. 

Before  proceeding  with  the  evaluation  of  the  source  terms,  it  will  be  convenient 
to  write  £[  as  the  product  of  a  frequency  term  and  a  term  that  varies  only  with  space 
and  time: 

£ i ( P . z , 0) , t )  =  q(p,z,t)  (3(to)  (25) 


where 


4(p.z,t)  = 


Anc(p,z,t) 

<nc) 


(26) 


is  a  random  function  of  the  electron  density  fluctuations  and 


p(co)  = 


<«n> 


Cl)  “  -  (  0)  2  '■ 


11) 


p/ 


is  a  deterministic  function  of  frequency  and  the  mean  free  electron  density. 
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For  the  first  source  term  in  Equation  23  a  straiehtforward  application  of  the 
Novikov  theorem  yields: 


S,  =  k i P ( co ! )  <^(p1,z,t1)  U(p , ,z,co, ,tj )  U*(p2,z,o>2,t2)) 


(28) 


OO  OO  OO 

=  k,p(0),)  Jdz'  J  d 2 p 1  Jdt'  <s(pi,z,ti)^(p',z’,t')> 

-  OO  -oo  -  oo 


X 


< 


SlKpi'Z.COj.t;) 


SU*(p2 


SflPWTT  lj’*<P2.Z.O)2,t2)  +  U(P,,Z.(0„I,) 


,z,co2,t2)\ 

SZ',0  /  • 


At  this  point  the  electron  density  fluctuations  are  assumed  to  be  stationary  and 
delta  correlated  along  the  z  axis.  This  Markov  assumption  has  the  mathematical  form 

(q(p,z,t)q(p',z',f)>  =  8(z-z')  A*(p-p',t-f)  (29) 

where  8()  is  the  Dirac  delta  function.  The  structure  function  A^(p-p',t-t’)  is  the 
autocorrelation  function  of  the  electron  density  fluctuations.  The  Markov  assumption 
is  discussed  in  some  detail  by  Tatarskii  (1971,  §64]  and  is  based  on  the  fact  that 
fluctuations  in  the  dielectric  constant  in  the  direction  of  propagation  have  little  effect 
on  the  transverse  fluctuation  characteristics  of  the  electric  Field.  It  is  the  fluctuations 
of  the  dielectric  constant  transverse  to  the  direction  of  propagation  that  dominate  the 
scattering  and  the  transverse  fluctuations  of  the  electric  field. 

Under  the  assumption  of  small-angle  scattering  for  which  the  parabolic  wave 
equation  is  valid,  the  component  of  the  electric  field  traveling  in  the  backward 
direction  will  be  negligible  compared  to  the  component  traveling  in  the  forward 
direction.  The  electric  field  U(p,z,co,t)  may  then  be  assumed  to  depend  on  ^(p',z',f) 
only  for  z’  <  z  (i.e.  the  electric  field  does  not  depend  on  electron  density  irregularities 
that  have  not  yet  been  encountered  along  the  forward  propagation  path).  Also, 
U(p.z,o),t)  depends  on  ^(p',z',f)  only  for  t'  <  t  (i.e.  the  electric  field  does  not  depend 
on  irregularities  that  have  not  yet  occurred).  Thus  SU(p,z,C0,t)/6^(p,,z,,t')  is  equal  to 
0  for  z'  >  z  and  for  t'  >  t. 

ddie  source  function  S]  may  now  be  rewritten  as 
z  °°  t 

Si  =  k,P,  J  dz'  6(z-z')  J  d 2 p ’  Jdf  A^(prp',trf)  x  (30) 

-OO  -oo  -OO 


/8U(pt,z,fa)t,tt) N  SU*(p2,z,o)2,t2)  \ 
\  8^(p',z',f)  U  (P2.z,«2,t2)  +  U(p1,z,co„t1)  / 


where  (3j  is  given  by  Equation  27  evaluated  at  (Oj.  Recalling  that 
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(31) 


x 

J  8(x’-x)  dx' 

-CO 


1 

1  ’ 


the  source  term  is  further  reduced  to 
»  t 


S,  =  ‘  J  d 2p’  Jdf  A^prp’.tj-t')  x 


(32) 


8U(pi,z,a)1,t1) 

8q(p',z,t') 


U*(p2,Z,G)2,t2)  + 


U(p(,z,(Oi,t,) 


5U*(p2,Z,CD2,t2)\ 

8^(p',z,t')  / 


The  parabolic  wave  equation  is  used  to  evaluate  the  functional  derivatives  6U/5^. 
Integrating  this  equation  from  -c*>  to  z  results  in 

z 

j  V  2  U(p,z",G),t)  dz"  -  2i(k)  (U(p,z,co,t)  -  u0(p,(0)]  -  (33) 

-oo 


Z 

<k>2(Ko)  J  ^(p>z",t)  U(p,z",co,t)  dz"  =  0 

-oo 


w'here  U0(p,m)  is  the  transmitted  signal.  After  applying  the  operator  5/5^(p’,z',t’), 
where  -°°  <  z  <  z  and  -°°  <  t'  <  t,  and  noting  that 

=  5<z-z'>  5<p-p')  sc-1')  ■  <34> 

Equation  33  becomes 

2i(k)  +  (k>2  P(co)  S(p-p')  S(t-f)  U(p,z,(0,t) 


+ 


z 

|  [  (k)2P(o))^(p,z",t)  -  vl 
z 


]  5U(p,z",co,t) 

J  S^(p',z’,t') 


dz” 


0  . 


(35) 


The  lower  limit  of  the  integral  in  this  equation  is  z'  because  8U(p,z,(o,t)/5£,(p',z’,t')  is 
zero  for  z  <  z .  Because  the  source  term  contains  factors  of  the  form 
5U(p,z,0),t)/8q(p',z,f ),  z  is  set  equal  to  z  in  Equation  35  with  the  result 


5U(p,z,0),t)  _  i(k)3(co) 
5q(p\z,t')  “  2 


5(p-p')  5(t-t')  U(p,z,co,t)  . 


(36) 
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After  substituting  this  into  Equation  32,  the  source  term  becomes 

OO  CO 

St  =  J  d2p'  J dt’  A^prp'.ti-O  r 

-OO  - OO 


ik i pi  _  ik^BiP? 

_-4~S(Pi-P')  S(ti-t')  -  -^L5(p2-p')  S(t2-f) 


(37) 


'Mi 


a§(0,0)  r  - A=(prp2,trt2)  r  . 


A  similar  expression  may  be  written  down  for  the  second  source  term  in  Equation  23: 

S2  =  k2p(co2)  <£(p2,z,t2)  Ufpj.z.coj.t!)  U*(p2,z,o)2>t2)>  (38) 

ik,(3ip2  ik^p2 

=  4  A^(prp2,trt2)  r  -  -f2  a^(0,0)  r  . 

2.2.4  Second  Form  of  Transport  Equation. 

The  transport  equation  for  the  mutual  coherence  function  is  now  given  by 
combining  Equations  37  and  38  with  23  with  the  result 


3T  i 

dz  +  2 


-i_T7  ^  —\7  ^ 

kj11  *  k2  ^ 


r 


(39) 


-  I  [2kik2p,p2  A^(p,-p2,t,-t2)  -  (kfp?  +  k2p2)  A0]  r  =  0 


where  A0  =  A^(0,0). 

The  differential  equation  for  T  will  be  solved  by  first  letting  T  equal  r0Ti  where 
T0  is  the  free  space  solution  to  the  transport  equation.  The  well-known  solution  of  the 
wave  equation  for  the  electric  field  in  free  space  (Eqn.  18  with  e,  =  0)  may  be  written 
down  directly.  The  Fresnel  approximation  that  z  »  Ipl  is  then  used  to  expand  the 
electric  field  and  the  free  space  solution  T0  is  computed.  The  quantity  T0  contains  the 

1/z2  term  that  partly  determines  the  mean  power  at  the  receiver.  The  next  step  is  to 
derive  a  differential  equation  for  Tj  from  Equation  39  and  the  free  space  solution.  It 
is  the  mutual  coherence  function  T]  that  determines  the  second-order  statistics  of  the 
received  signal. 


2.2.4. 1  Free  Space  Solution  Tn  In  free  space  and  for  spherical  geometry,  the 
complex  envelope  of  the  electric  field  is 


E  =  E, 


exp  f-iklrl) 


(40) 
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It  is  easy  to  verify  that  this  is  a  solution  of  the  wave  equation  with  c,  set  to  zero. 
Under  the  assumption  of  small-angle  scattering,  z2  »  x2  +  y2  in  the  region  of  interest 
and  Irl  may  be  expanded  as 

_ _ 2  2 

/  s  s  x  +  y 

Irl  =  v  x“  +  y“  +  zY  =  z  + — —  .  (41) 

With  this  Fresnel  approximation, 

-  Eo  r  . ,  r  .  X2  +  y2  \] 


E  =  yexpj-  ik  |z  +  jj  .  (42) 

After  recalling  that  U  =  exp  (ikz)  F,  the  free  space  mutual  coherence  function  is 

r„  =  < U ( p , , /.,(!),)  U*(p2.z,(o2)>  (43 ) 

1  |  ik,(x]  +  Vj)  ikdx4  +  v2)  "] 

=  :.2exP‘  ;  . 


2. 2. 4. 2  Differential  Equation  for  I',.  After  substituting  V  =  1 3,1' s  into  the 
•ort  equation  and  using  Equation  43  for  the  free  space  solution,  the  equation  for 


T,  becomes 

dr i  j  ->  i  2 

dz+zk^1  (44) 

.  i  »i  +  y±  <K±  +  *4  *1± .  y 2  »1±  c  r  _  () 

[_  z  3x!  z  dyt  z  dx2  z  dy2j  '  1 

where  the  source  term  is 

^  =  S  t“k  tk’PtPz  ^sfPi-P-'ti't’)  -  (ktPt  +  k  2P2)  A(>  ]  (45) 

In  order  for  to  represent  a  statistically  stationary  random  process  in  space, 
frequency,  and  time,  T,  must  be  a  function  only  of  the  differences  p , -p and 
t| -t2  .  It  is  therefore  useful  to  transform  liquation  44  to  sum  and  difference  spatial  and 
frequency  coordinates: 


v;  = 


X|  +  X2 

Y 

Vi  +  y 

=  -> 

X|  -  X2 

n  : 

=  yi  -  y: 

r)2  d2 

■> 

d2 

dX2  +  dY2 

vd 

-  -^2  + 

VS'Vd 


d2  d2 

clXclC  +  bYr>r| 


ks 


k  L +Jk  ’ 


k(j  —  k |  "  k2 


After  some  manipulations,  the  equation  for  fj  reduces  to 


^r, 

(h 


4k^  -  k  d 


kd^d  +  4  Vs  ■  ksVs  Vd  '  T]  + 


(47) 


1 

/ 


+ 


-sr, 


=  o 


The  boundary  condition  for  this  equation  is  that  Tj  evaluated  at  z  =  -zt  must  be 
equal  to  unity  independent  of  the  other  spatial  coordinates.  Also,  the  source  term  S 
under  most  conditions  is  a  function  only  of  the  difference  coordinates.  It  will 
therefore  be  assumed  that  T,  is  independent  of  X  and  Y.  However,  the  source  term 
will  be  a  function  of  X  and  Y  if  the  spatial  extent  of  the  scattering  region  is  small  as, 
for  example,  in  a  barium  cloud.  The  assumption  that  T,  is  independent  of  X  and  Y 
Lhen  requ'res  that  the  disturbed  region  in  the  ionosphere  be  large  compared  to  the 
region  from  which  scattered  signal  energy  is  received. 

The  transport  equation  may  be  further  reduced  by  noting  that  the  1/z  term,  when 
/  is  large,  will  be  small  compared  to  the  other  terms  and  may  be  neglected.  'The 
transport  equation  for  F]  then  becomes 


r)T, 

dz 


:>kd 

4k;-ka 


Vu  F,  -  ST,  =  0 


(48) 


2.3  DKI/I A-LAYKR  APPROXIMATION. 


As  an  RF  wave  propagates  through  a  thick,  irregularly  structured  ionization 
layer,  the  wave  first  suffers  random  phase  perturbations  due  to  random  variations  in 
the  index  of  refraction.  As  the  wave  propagates  farther,  diffractive  effects  introduce 
fluctuations  in  amplitude  as  well  as  phase.  If  the  standard  deviation  of  the  phase 
fluctuations  that  are  suffered  by  the  wave  is  large,  then  the  amplitude  fluctuations  are 
characterized  by  a  Rayleigh  probability  distribution  when  the  wave  emerges  from  the 
layer.  The  delta  layer  approximation  assumes  that  die  phase  and  amplitude  fluctuations 
are  imparted  on  the  wave  in  an  infinitesimally  thin  layer.  This  assumption  is  consistent 
with  the  Maikov  assumption  made  in  deriving  the  differential  equation  for  the  mutual 
coherence  function. 


An  analytic  solution  for  the  two-position,  two-frequency  mutual  coherence 
function  has  been  obtained  for  plane  waves  by  Sreenivasiah,  Ishimaru,  and  Hung 
|197bi  for  a  thick  ionization  layer.  Wittwcr  |]97A)  extended  this  solution  to  tr  at 
spherical  waves.  The  analytic  form  of  this  solution  is  sufficiently  complex,  however, 
that  the  necessary  Fourier  transforms  required  to  compute  the  CiPSD  cannot  be 
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performed  in  closed  form.  The  complex  analytic  form  is  simplified  by  the  use  of  the 
delta-layer  approximation  to  obtain  tractable  expressions  for  the  mutual  coherence 
function  and  the  GPSD.  Wittwer  has  evaluated  the  accuracy  of  the  delta-layer 
approximation  as  it  affects  the  delay  distribution  of  the  received  signal  and  has  found 
that  the  maximum  error  is  small  for  transionospheric  satellite  communication  link 
geometries  as  long  as  the  parameters  of  the  GPSD  are  properly  selected.  Wittwer 
[19791  has  derived  expressions  for  the  signal  parameters  of  the  GPSD  that  include  the 
effects  of  a  thick  scattering  layer. 

Now  a  relationship  between  the  electron  density  fluctuations  and  the  phase 
variations  imposed  on  the  wave  may  be  calculated.  The  differential  phase  change  of 
the  wave  along  the  propagation  path  l  is 

=  rc  X  Anc(p,z,t)  .  (49) 

Under  the  assumption  of  small-angle  scattering  used  to  derive  the  differential  equation 
for  T | ,  d/di  is  approximately  equal  to  d/dz,  and  the  total  phase  change  of  the  wave  is 


4>  =  rcX  J  Anc(p,z,t)  dz  =  rcU(nc)  j  £(p,z,t)  dz  ,  (50) 

integrated  through  the  ionization  layer.  The  autocorrelation  function  of  the  phase 
change  is 

<<J)(p,t)<t>(p\t')>  =  [re>w<ne)]2  JdzJdz'  <£(p,z,t)£(p\z\t)> 


=  [rcA<nc>]2  J  A^(p-p\t-t')  dz 


(51) 


=  [rcA.(ne>]2  L§  A^(p-p',t-t’) 

where  L§  is  the  delta  layer  thickness.  The  Markov  approximation  (Eqn.  29)  has  been 
used  in  evaluating  the  autocorrelation  of  However,  it  is  show-n  m  Appendix  A  that 
Equation  51  is  valid  as  long  as  the  scattering  layer  thickness  is  large  compared  to  the 
parallel  decorrelation  distance  of  the  electron  density  fluctuations.  The  phase  variance 
imj  urted  on  the  wave  is 

=  [re^(nc>]2  U  A0  .  (52) 

The  quantity  A0  depends  on  the  power  spectrum  of  the  electron  density  fluctuations  in 
the  ionosphere.  The  value  of  A„  for  a  three-dimensional  K  4  in  situ  power  spectrum 
and  for  the  delta-layer  approximation  is  given  in  Appendix  B. 
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In  general,  only  part  of  the  total  phase  variance  results  in  the  Rayleigh  amplitude 
fluctuations  that  are  "described  by  the  GPSD.  Wittwer  [1979,  1980]  calls  this  part  the 
Rayleigh  phase  variance.  The  rest  of  the  total  phase  variance  is  associated  with  the 
mean  dispersive  effects  described  by  the  exponential  term  of  Equation  17.  It  is  the 
smaller-sized  electron  density  fluctuations  that  result  in  diffractive  effects  and  the 
larger-sized  fluctuations  that  result  in  dispersive  effects.  Wittwer  [1982]  describes  how 
the  Rayleigh  phase  variance  may  be  separated  from  the  total  phase  variance.  In  the 
developments  that  follow,  the  phase  variance  in  Equation  52  will  be  assumed  to  be  the 
Rayleigh  phase  variance  associated  with  diffractive  effects  and  amplitude  fluctuations. 

2.4  FORMAL  SOLUTION  OF  THE  TRANSPORT  EQUATION. 


Before  proceeding  with  the  solution  of  the  transport  equation,  it  is  convenient  to 
expand  the  source  term  S  by  making  two  non-re strictive  assumptions.  First,  it  will  be 
assumed  that  the  RF  frequencies  of  interest  are  large  compared  to  th*  plasma 
frequency.  Second,  it  will  be  assumed  that  k<i  is  much  smaller  than  ks.  The  solution 
obtained  will  then  be  valid  for  a  small  range  of  frequencies  around  ks  and  for 
frequencies  large  compared  to  typical  peak  ionospheric  plasma  frequencies  of  a  few 
hundred  MHz.  With  these  assumptions,  the  source  term  becomes 


S  = 


(COp) 

4c4k 


2 

I  [A^(p,t)  -  A0]  - 
s 


kd(tQp)2 
8c4k4  A° 


(53) 


where  it  is  understood  that  p  is  a  two-dimensional  relative  position  vector  in  the  ^-r| 
plane  normal  to  the  line-of-sight,  and  t  is  time  difference. 


The  last  term  in  the  equation  for  S  is  a  function  of  frequency  and  z  but  is 
independent  of  ^  and  r\.  This  suggests  that  another  useful  factorization  is  Tj  =  r2r3 
where  T3  is  independent  of  £  and  r\.  After  making  this  substitution  in  the  transport 
equation  for  T[  and  separating  variables,  the  result  is 


1  ar2 
T2  dz 


ikd  Vfo 
2k  I  r2 


(<Dn>2  1  dr,  kj<0)p>2 

4A|  ia^(p'0  ‘ Ao)  +  r,  3z  +  8?i4  A° 


0  .  (54) 


The  last  two  terms  of  this  equation  depend  only  on  the  spatial  variable  z.  Therefore 
the  sum  of  these  two  terms  must  be  separately  equal  to  zero  for  arbitrary  values  of  the 
other  spatial  variables  C,  and  T|. 


The  source  term  in  the  T3  equation  is  only  non-zero  within  the  delta  layer.  Thus 
from  the  transmitter  to  the  delta  layer,  T3  is  unity.  With  this  boundary  condition,  the 
solution  of  the  f3  equation  is 


T3  =  exp 


kd<0)P>2Ao  ,7  7  . 
8c4k,  (  0 


,  z  >  zt  . 


(55) 


This  term  gives  the  effect  of  different  transit  times  of  different  frequencies  that  results 
from  the  frequency  dependence  of  the  index  of  refraction. 
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Now  the  equation  for  T2  may  be  solved  using  the  delta-layer  approximation. 
The  equation  for  f2  is 


1  dr2  ikd  VdT2  ((Ofr)" 
T2  dz  2ks  T2  4cJk; 


|A*(p,t)  -  A0) 


=  0  . 


(56) 


Within  the  delta  layer,  the  k^  term  is  small  compared  to  the  source  term  and  may  be 
ignored.  This  is  equivalent  to  assuming  that  diffractive  effects  are  not  important 
within  the  delta  layer  [Lee  and  Jokipii  1975b],  Integrating  Equation  56  through  the 
delta  layer  gives  the  value  of  r2  at  the  point  where  the  wave  emerges  from  the  delta 
layer: 


j  («2n>2  1 

r2  =  exp  |  4cPk2  [A'(p.t)  -  A0]  [ 


(57) 


The  solution  from  this  point  proceeds  as  follows.  Between  the  delta  layer  and 
the  receiver,  the  signal  propagates  in  free  space  so  the  equation  for  T ,  is 


dr 


i  _2ikd  y 2  r 
dz  '  4kg-kj  Vd  r‘ 


0  . 


(58) 


“> 

The  V^T,  term  in  this  equation  gives  the  effects  of  diffraction  on  the  signal  as  it 
propagates  from  the  delta  layer  to  the  receiver.  The  boundary  condition  is  T,  = 
at  the  point  where  the  wave  emerges  from  the  delta  layer.  This  equation  is  easily 
solved  by  taking  the  Fourier  transform  from  spatial  coordinates  C,  and  r|  to  angular 
coordinates.  First,  it  is  convenient  to  transform  variables  to 


u 


and  v 


n 

z  ‘ 


(59) 


The  angular  variables  Ku  and  Kv  are  then  independent  of  z.  After  the  change  in 
variables,  the  equation  for  T,  becomes 


dr,  ikd  Id2  d21_  n 
dz  -  2k  2z2  Ldu2  +  dv2J  V]  -  ° 


(60) 


Fourier  transform  pairs  from  spatial  coordinates  p  in  the  plane  normal  to  the 
line-of-sight  to  angular  coordinates  Kj_,  from  carrier  frequency  differences  0)  to  delay 
x,  and  from  time  differences  t  to  Doppler  frequencies  md  are  defined  in  this  report  to 
be 
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(61a) 


r(K±)  =  j  exp(-iKi  p)  Hp)  d2p 


A  r  dW 

r(T)  =  J  exp(icot)  T(co)  ^ 


(61b) 


r ( co D )  =  J  exp(icoot)  r (t)  dt 


(61c) 


r  a  d‘Ki 

Hp)  =  J  exp(iKi-p)  T(Ki)  ^2 


(62a) 


r(w)  = 


J  exp(-i 


icox )  f(T)  dt 


(62b) 


f  a  d(i)n 

r(t)  =  J  exp(-KODt)  r(WD)  2K 


(62c) 


Upon  transforming  from  u  and  v  to  Ku  and  Kv,  the  equation  for  f  i  is 

A 

^-^(Ki+Ki)r.  -  0  . 


This  equation  is  integrated  from  z  =  zt  +  L5  to  z  =  z\  +  zr  which  gives 


a  a  1Ndr/  7  7  \ 

r,(Ku,Kv,Zt+zr,co,t)  =  ri(Ku,Kv,zt+L5,co,t)  exp^--^y(Ku  +  Kvj 


where 


Zt+Zr 

fdJ 


rC 

J 

Zt+Ls 


Zr  -  L5 
(Zt+Zr)(Z(+L§) 


The  expression  for  y  may  be  simplified  by  setting  Lg  to  zero,  which  is  consistent  with 
the  delta-layer  approximation.  Thus  7  is  given  by 


(66) 


_ Zr  _ 

^  ~  zt(zt  +  Zr)  ‘ 

At  this  point,  a  formal  solution  for  the  Fourier  transform  of  the  two-position, 
two-frequency,  two-time  mutual  coherence  function  at  the  receiver  has  been  computed 
in  terms  of  the  structure  function  A^(p,t)  -  A0.  In  the  next  two  sections,  models  for 
the  temporal  variation  of  the  structure  function  and  the  quadratic  approximation  for 
the  spatial  variation  of  this  function  will  be  discussed. 

2.5  TEMPORAL  VARIATION  MODELS  FOR  Ane. 

Under  the  frozen-in  model,  the  temporal  variation  of  the  electron  density 
fluctuations  is  given  by 

Anc(p,z,t)  =  Anc(p-vt,z,0)  .  (67) 

This  equation  is  valid  if  the  electron  density  fluctuations  with  a  scale  size  L0  do  not 
appreciably  change  their  shape  within  the  time  required  for  the  structures  to  drift  a 
distance  L0.  This  is  called  Taylor's  fiozen-in  hypothesis,  and  this  model  has  been  used 
in  most  previous  calculations  of  the  mutual  coherence  function  for  RF  propagation 
through  the  ionosphere.  Using  this  model,  the  structure  function  A  becomes 

A^(p,t)  =  A^(p-vt,0)  .  (68) 

This  model  is  accurate  for  ionospheric  conditions  where  the  ionization  has  broken  up 
into  a  single  layer  of  striations  aligned  with  the  earth's  magnetic  field  lines. 

The  frozen-in  model  may  not  be  valid  before  striations  have  formed  or  when 
there  are  multiple  scattering  layers  in  the  line-of-sight.  Under  these  conditions,  a 
“turbulent  model”  that  decorrelates  spatial  and  temporal  variations  may  be  more 
appropriate. 

Dr.  Leon  A.  Wittwer  of  the  Defense  Nuclear  Agency  has  recently  proposed  a 
more  “General  Model”  for  the  temporal  and  spatial  behavior  of  the  structure  function 
[Wittwer  1988].  In  this  model  it  is  assumed  that  the  stnm.ure  function  depends  on  time 
and  space  as  p  -  Cvt,  where  C  is  the  space-time  correlation  coefficient  of  the  electron 
density  fluctuations.  This  model  then  varies  smoothly  between  the  frozen-in  model 
where  C  is  unity  and  the  turbulent  model  where  C  is  zero. 

2.6  QUADRATIC  PHASE  STRUCTURE  APPROXIMATION. 

As  a  practical  matter,  the  quadratic  phase  structure  approximation  is  required  to 
make  the  exponent  in  Equation  57  quadratic  in  the  spatial  and  temporal  variables.  This 
gives  the  resulting  mutual  coherence  function  a  tractable  mathematical  form. 
However,  for  the  small  angle,  strong-scattering  conditions  considered  in  this  report, 
the  correlation  distance  of  the  signal  will  be  much  smaller  than  the  correlation  distance 
of  the  electron  density  fluctuations.  The  mutual  coherence  function  will  then  be 
determined  primarily  by  the  values  of  A^(p,t)  -  A0  at  small  values  for  the  space  and 
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time  differences,  and  a  Taylor  series  expansion  of  A^(p,t)  -  A0,  keeping  only  the 
quadratic  terms,  will  provide  a  reasonably  accurate  mutual  coherence  function. 

First  consider  only  the  spatial  dependence  of  the  structure  function,  A^(p).  The 
temporal  behavior  will  be  included  later.  The  Taylor  series  expansion  will  make  the 
following  calculations  independent  of  the  functional  form  of  A^(p)  as  long  as  the 
'£cor.d  derivative  of  A^(p)  exists  for  p  equal  to  zero.  A  detailed  discussion  of  A*(p) 
will  therefore  be  deferred  to  Appendix  B.  However,  some  generic  properties  of  A^(p) 
need  to  be  considered  here  in  order  to  specify  the  functional  form  of  the  quadratic 
terms  in  the  Taylor  series  expansion. 

The  propagation  coordinate  systems  are  shown  in  Figure  2.  The  z  axis  is  along 
the  line-of-sight  and  the  t  axis  is  aligned  with  the  geomagnetic  field  lines  B  at  the 
elevation  of  the  delta  layer.  The  r  and  s  axes  are  in  a  plane  normal  to  the  magnetic 
field.  The  penetration  angle  O  is  the  angle  between  the  z  axis  and  the  t  axis.  At  the 
receiver,  the  x  axis  direction  is  given  by  the  cross  product  of  B  and  the  line-of-sight 
unit  vector,  z,  and  the  y  direction  is  orthogonal  to  both  x  and  z. 

Now  consider  the  functional  form  the  power  spectrum  S^(K)  of  the  electron 
density  fluctuations  in  the  r-s-t  coordinate  system.  The  power  spectrum  is  usually 
assumed  to  be  a  function  of  the  quantity 

Lr2Kr2  +  L2K2  + 


sin  O 


Figure  2.  Propagation  coordinate  systems. 
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where  Lr ,  Ls ,  and  Lt  are  the  scale  sizes  of  the  fluctuations  in  the  three  directions.  The 
usual  assumption  that  the  electron  density  fluctuations  are  elongated  in  the  t  direction 
and  are  symmetric  about  the  t  direction  will  be  made.  The  scale  sizes  are  then 

Lr  =  Ls  =  L0  (69) 

Lt  =  q  Ln 

where  q  (q  >  1 )  is  the  axial  ratio. 


The  power  spectrum  S*(K)  will  now  be  evaluated  in  the  x-y-z  coordinate  system 
in  order  to  calculate  A=(p)  in  this  system.  The  power  spectrum  in  x-y-z  coordinates  is 
a  function  of  the  quantity  [Wittwer  1979) 


L2K2  +  L~yKl  +  L;K2a  +  2LyzKyKz 


where 

Lx  =  L() 

Ly  =  L0  Vcos2d>  +  q2sin20  (70) 

Lz  =  L0  V  sin20  +  q2cos2<t> 

Ly z  =  L()  (q2-l)  sinO  cosO 

The  structure  function  A*(p)  is  the  two-dimension  Fourier  transform  of  the 
power  spectrum  S=(K): 


A*(p)  =  j  exp  (iKi  p)  S^(Ki,Kz=0) 


d2Ki 

{2k)2 


(71 


Setting  Kz  =  0  in  this  equation  results  in  S^(K)  being  a  function  of 


L2K2  +  LyK2 


and  the  structure  function  being  calculated  in  the  plane  of  the  delta  layer  (i.e.  at  z  =  0). 
After  the  Fourier  transform  is  performed,  A=(p)  will  be  a  function  of  the  quantity 


The  quadratic  phase  structure  approximation  then  takes  the  form 


A  =  (p)  =  A„  1 


A  - 


L: 


+ 


Ly 


\  ' 


I 


ill) 
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The  coefficient  A2  is  calculated  in  Appendix  B  for  a  K  4  electron  density  fluctuation 
spectrum.  The  coefficient  A0  is  given  in  terms  of  the  phase  variance  by  Equation  52. 


Under  the  general  model,  the  temporal  variation  of  the  structure  function  is  due 
to  a  translation  of  the  spatial  coordinates.  Equation  72  then  becomes 


A^(p,t)  -  A0 


1  -  Ai 


xi  yl  r_  Cxt  xt  Cvt  y  t 
,  2  +  ,  2  +  T2  '  ^  -l 


1  0 


LxTo 


Ly  Tq 


(73) 


The  quantity  T0  can  be  thought  of  as  a  “decorrelation  time”  of  the  electron  density 
fluctuations  in  the  ionosphere.  The  degree  of  space-time  correlation  is  determined  by 
the  coefficients  Cxt  and  Cyt . 

2.7  MUTUAL  COHERENCE  FUNCTION. 


The  boundary  condition  for  the  mutual  coherence  function  T,  at  the  delta  layer 
is  given  by  the  product  r2T3.  Using  Equation  52  to  relate  A,,  to  the  phase  variance  and 
then  setting  the  delta  layer  thickness  L§  to  zero  gives  the  following  for  V2  and  f3  at  the 
point  where  the  wave  emerges  from  the  delta  layer: 


r,  = 


exp 


-°JA2  *2 

l^x 


r 


y.2. 

+  L2  +  T2 

Ly  1  0 


Cxtxt  c^yt 

LxT0  LyTo 


and 


r,  = 


exp 


1  fG^COV' 

«o  J  J 


(74) 


(75) 


where  (o0  =  cks  is  the  radian  frequency  of  the  carrier. 

Upon  changing  variables  in  the  equation  for  T2  to  the  dimensionless  spatial 
cooiuinates  u  and  v,  performing  the  Fourier  transform  to  Ku  and  Kv  coordinates,  and 

A 

assuming  that  the  delta  layer  thickness  is  zero,  the  boundary  condition  for  T ,  is 


T  i(Ku,Kv,Zt,co,t) 


x  exp 


TtUxLy 

-  c4a2z2  ^ 

-  <4a2(cx2(  + 


2  .  II 

-  a4)A2  2 

1  o 


2  ...2  .  2  vl  ,  I  2y2 


2C0n 


LJM  +  LyK$ 

4o|A2Zi 


n  2  \  t2  .  CXlLXKU  +  CytLyKy  Jt 
Lyt/ 


rrj'2  ^ 

*0 


Tn 


(76) 


The  solution  at  the  plane  of  the  receiver  is  given  by  Equation  64  and  this  boundary 
condition. 
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After  performing  the  inverse  Fourier  transform  on  the  solution  T^Ku.Kv.Zt.w.t) 
and  converting  to  unnormalized  distance  units  x  and  y,  the  two-position,  two- 
frequency,  two-time  mutual  coherence  function  is 


T^x.y.toT)  = 


1 

hVl 

exp 
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The  decorrelation  distances  ix  and  !y  are  given  by  the  expressions: 
(zt  +  zr)Lx 

Jtv  = - -  j— : 

'  A 
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7  z^VXT 


zta^VA^ 

the  decorrelation  time  is  given  by  the  expression 

^0  j  i 

o<t>VA2 

and  the  coherence  bandwidth  (Oc0h  is  given  by  the  expression 

AyCOpL2  zt  +  zr 
(°coh  2ccj|a2  Ztzr  • 


(77) 
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Under  the  delta-layer  approximation,  the  ratio  of  the  two  decorrelation  distances 
is  a  function  of  only  the  penetration  angle  and  the  axial  ratio  and  is  given  by 


=  r*  =  V cos2<P  +  q2sin2$  . 

*X  LX 


(82) 


The  orientation  of  the  x  and  y  axes  is  usually  chosen  so  that  ix  is  the  smaller  of  the  two 
decorrelation  distances.  The  minimum  decorrelation  distance  is  also  referred  to  as  t0 
in  the  literature.  This  notation  will  occasionally  be  used  in  this  report.  Equation  82  is 
used  to  calculate  ly  when  only  the  minimum  decorrelation  distance  is  specified. 
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It  will  be  seen  in  the  next  section  that  the  coherence  bandwidth  is  proportional  to 
2k{q  where  f0  is  the  frequency  selective  bandwidth  of  the  signal.  It  is  clear  from  the 
form  of  the  mutual  coherence  function  in  Equation  77  that  the  coherence  bandwidth 
could  have  also  been  defined  by  either  of  the  products  coCohAy  or  coCohAx.  The 
asymmetry  factors  Ax  and  Ay  have  been  included  in  the  definition  of  the  coherence 
bandwidth  in  order  to  simplify  the  relationship  between  0)Coh  and  f0. 

The  equations  for  the  decorrelation  time,  decorrelation  distance,  and  the 
coherence  bandwidth  are  only  valid  under  the  delta-layer  approximation  and  do  not 
reflect  how  these  parameters  are  actually  calculated  for  a  given  ionospheric 
disturbance.  Wittwer  [1979]  has  derived  expressions  for  the  decorrelation  distance  and 
the  coherence  bandwidth  that  are  valid  for  more  general  scattering  layer  geometries. 
It  is  these  more  general  expressions  that  are  used  in  signal  specifications  to  calculate 
these  signal  parameters. 

For  satellite  communications  links,  the  distance  from  the  ionosphere  to  the 
satellite  is  typically  much  larger  than  the  distance  from  the  ionosphere  to  the  ground. 
If  the  transmitter  is  on  the  satellite  then  zt  »  zr  and  if  the  transmitter  is  on  the  ground 
then  zt  «  zr  .  Because  the  expressions  for  the  decorrelation  distances  are  not 
reciprocal  in  zt  and  zr,  the  values  of  £x  and  4y  depend  on  the  direction  of  propagation. 
However,  the  expression  for  the  coherence  bandwidth  is  reciprocal  in  zt  and  zr  so  (0coh 
and  the  frequency  selective  bandwidth  are  independent  of  the  propagation  direction.  It 
is  interesting  to  note  that  the  expression  for  the  decorrelation  time  is  independent  of  the 
propagation  distance,  depending  only  on  the  “decorrelation  time”  of  the  striations  and 
the  standard  deviation  of  the  phase  fluctuations. 

The  decorrelation  distance  of  the  signal  as  it  emerges  from  the  delta  layer  is 
approximately  equal  to  Lo/a^.  Under  strong-scattering  conditions  where  »  1  ,  this 
distance  will  be  much  smaller  than  the  electron  density  fluctuation  scale  size  and  the 
quadratic  phase  structure  approximation  conditions  are  met.  Conversely,  under  weak 
scattering  conditions  where  Cty  =  1,  the  quadratic  phase  approximation  will  give 
inaccurate  results  for  the  mutual  coherence  function. 

Now  consider  the  mutual  coherence  function  for  two  positions,  one  frequency 
and  two  times: 


T  i(x,y,0,t) 


n2 

ToJ. 


(83) 

The  quantities  ix,S.y,  and  i()  are  seen  to  be  the  1/e  points  of  T^x, 0,0,0), 
r^CfyAO),  and  T^O.O.O.t)  respectively  and  are  therefore  consistent  with  the  usual 
definitions  of  decorrelation  distance  and  time. 
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2.8  GENERALIZED  POWER  SPECTRAL  DENSITY. 

The  generalized  power  spectral  density  S(Kj_,t,0l>d)  of  the  signal  incident  on  the 
plane  of  the  receiver  is  the  Fourier  transform  of  the  mutual  coherence  function: 

oo  oo  oo 

S(Ki,t,cl>d)  =  /d2p  I  2k  ldt  rKP-w-0  exp  [-i(Ki  p-coT-WDt)]  (84) 

-OO  -OO  -oo 

where  angle-of-arrival  Ki  is  the  Fourier  transform  pair  of  position  p  in  the  x-y  plane, 
delay  t  is  the  Fourier  transform  pair  of  relative  carrier  frequency  to,  and  Doppler 
frequency  cod  is  the  Fourier  transform  pair  of  relative  time  t.  The  quantity 
S(Ki,T,C0D)(d2Kiy4jr2)dT(dc0D/2rt)  is  equal  to  the  mean  signal  power  arriving  with 
angles-of-arrival  in  the  interval  Kj_/4tc2  to  (K_L+d2Ki)/4rc2,  w-ith  delays  relative  to  a 
nominal  propagation  time  in  the  interval  t  to  i+di,  and  with  Doppler  frequencies  in  the 
interval  coo/2Jt  to  (cod  +  dcOD)/2Jt.  The  delay  dependence  of  the  GPSD  is  a 
consequence  of  the  fact  that  some  of  the  signal  energy  takes  a  dog-leg  path  through  the 
ionosphere  from  the  transmitter  to  the  receiver  and  arrives  later  than  the  signal  energy 
that  propagates  straight  through  the  ionosphere. 

In  general,  the  GPSD  can  be  written  as  the  product  of  a  Doppler  spectrum 
Sd(c*)d)  and  an  angle-delay  spectrum  Skt(K±,t): 

S(Ki,t,cod)  =  Sd(«d)  Skt(Ki,t)  .  (85) 

After  performing  the  integrals  indicated  in  Equation  84,  the  Doppler  spectrum  for  the 
general  model  is 


(T0WD  -  C*tKX  '  Cy(K  v &v  )“ 

40  •  C,\  -  C),)  ,  ■  <86) 

(General  Model) 

The  angle-delay  part  of  the  GPSD  for  the  general  model  is  identical  to  that 
derived  previously  for  the  frozen-in  and  turbulent  models  [see,  for  example,  Dana 
1986].  In  terms  of  the  components  of  Kj_  (Kx  and  Ky)  the  angle-delay  part  of  the 
GPSD  is 
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where  the  delay  parameter  a  is  defined  to  be 


(87) 


24 


a 


G)(, 

Cy<j>G\oh 


(88) 


The  components  of  are  related  to  the  scattering  angles  0X  and  0y  about  the  x  and  y 
axis  respectively  by  the  relations 


^  2k  sin(0x) 
Kx  -  ^ 


(89a) 


2k  sin(0y) 
K>  -  \ 


(89b) 


It  should  be  noted  that  the  range  of  delay  in  Equation  87  is  from  -«>  to  +<». 
However,  this  delay  is  relative  to  some  nominal  propagation  time,  and  the  value  of  the 
GPSD  rapidly  approaches  zero  for  decreasing  negative  values  of  G)Coht.  Thus  negative 
values  of  delay  present  no  real  problem  with  causality. 

It  is  interesting  to  examine  the  limits  of  the  general  model  Doppler  spectrum  in 
order  to  show-  that  this  model  does  indeed  encompass  both  the  frozen-in  and  fully 
turbulent  models.  These  limits  for  the  Doppler  spectra  are: 

Limit 

Cxt-  1  Sq(wd)  =  2k  t0  8(t0cod  -  Kx£x)  (Frozen-in  Model)  (90) 

Cy{  — *  0 


and 


Limit  r  ^^21 

CXt-*0  Sd(md)  =  V^T()exp  — !  •  (Turbulent  Model)  (91) 

C  y  t  — »  0  L  4  j 

For  the  frozen-in  model,  the  delta-function  relationship  between  Doppler  frequency 
and  Kx  is  what  is  obtained  by  assuming  that  the  random  diffraction  pattern  of  the  signal 
is  “frozen”  and  drifts  in  the  x-direction  past  the  receiver.  For  the  turbulent  model,  the 
Doppler  spectrum  is  independent  of  Kx  so  the  temporal  and  spatial  variations  in  the 
received  signal  are  also  independent.  This  Gaussian  form  is  the  true  turbulent  Doppler 
spectrum,  as  opposed  to  previous  turbulent  models  that  assumed  an  gjd  Doppler 
spectrum  for  mathematical  convenience  [Dana  1986]. 

The  significance  of  the  parameters  a,  £x,  iy,  and  (Ocoh  that  appear  in  the  angle- 
delay  part  of  the  GPSD  will  be  discussed  in  the  next  subsections. 

2.8.1  Delay  Spread  and  a. 

The  distribution  of  delayed  signal  power  with  a  fixed  angle-of-arrival,  which  is 
given  by  the  second  exponential  term  in  Equation  87,  has  the  Gaussian  form 
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Sx(T) 


acocoh  r  «2^coh  (T  -  tp)2 

4Tk  exp  L"  2 


(92) 


where  tp  is  the  additional  propagation  time  for  signals  arriving  at  the  angles  Kx  and 
Ky.  To  see  this,  the  expression  for  tp  is  expanded  using  the  definitions  for  4X,  4y,  and 
0>coh  and  the  assumption  of  small-angle  scattering  to  give 

_  9x  +  ey  zr(zt  +  zr) 
lP  "  2c  z,  • 

The  geometry  of  the  scattering  in  one  plane  is  shown  in  Figure  3.  Clearly,  for  small- 
angle  scattering,  the  angle  9t  at  the  transmitter  is  related  to  the  receiver  scattering  angle 
0r  by  the  relationship 


6t  =  7  Br  •  (94) 

zt 

The  difference  d  between  the  line-of-sight  distance  and  the  scattered  path  length  is 
given  by 

,  FI  TZ  n  TT  ,  ,  9r  zr(Zt  +  z r) 

d  =  V  ^  +  z  r0r  +  V  zt  +  z  t0[  -  (zr  +  Zt)  =  J  z 

for  scattering  in  one  plane  only.  When  scattering  about  both  the  x  and  y  axes  is  taken 
into  account,  the  total  path  difference  is  given  by  Equation  95  with  replaced  by 
0^-0y.  The  additional  time  required  for  the  signal  to  propagate  along  the  scattered 
path  is  just  d/c  which  is  equal  to  tp. 


Scattering 

Layer 


Figure  3.  Scattering  geometry. 
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For  a  given  vaiue  of  the  coherence  bandwidth,  larger  values  of  a  result  in  a 
narrower  delay  distribution  of  signal  power  about  the  time  tp.  The  delay  parameter  a 
is  then  a  measure  of  the  relative  importance  of  diffraction  and  dispersion  with  large 
values  of  a  indicating  strong-scattering  effects  and  small  values  indicating  weak 
scattering  or  dispersive  effects.  The  strong-scattering  limit  then  requires  that  the  value 
of  a  be  large. 

2.8.2  Frequency  Selective  Bandwidth  and  coCOh- 

The  frequency  selective  bandwidth  f0  is  an  important  measure  of  the  effects  of 
scintillation  on  the  propagation  of  wide  bandwidth  signals.  This  quantity  is  defined  as 
the  standard  deviation  of  the  time-of-arrival  jitter  ax: 


fo 
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27ta-c 


(96) 


where 

a?  =  (t2>  -  <t>2  .  (97) 

These  delay  moments  can  be  calculated  directly  from  the  angle-delay  part  of  the  GPSD 
using  the  equation 


Po  <xn> 
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(98) 


It  is  easy  to  show  that  the  mean  signal  power  P0  is  equal  to  unity.  The  first  and  second 
moments  are  easily  obtainable  in  closed  form  giving  the  relationship  between 
coherence  bandwidth  and  the  frequency  selective  bandwidth: 


t^coh  -  2ltfo 
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-i  1 
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(99) 


This  expression  is  valid  under  the  quadratic  phase  structure  approximation  for  which 
the  GPSD  is  valid.  Yeh  and  Liu  [  1 977}  have  calculated  an  expression  for  the  time 
delay  jitter  keeping  both  the  second  and  forth  order  terms  in  the  expansion  of  the 
structure  function  A^(p).  Their  result  has  more  terms  in  the  expression  for  the  time 
delay  jitter.  However,  the  additional  terms  are  significant  only  when  the  quadratic 
approximation  for  the  structure  function  is  invalid  and  therefore  only  when  the  GPSD 
derived  above  is  also  invalid. 

The  1  +  1  /a2  term  in  the  expression  for  the  coherence  bandwidth  represents  the 
relative  contributions  to  the  time  delay  jitter  from  diffraction  (1)  and  dispersion 
(1/a2).  In  the  limit  that  a  is  large,  the  time  delay  jitter  is  determined  by  diffractive 
effects  alone  which  should  be  the  case  under  strong-scattering  conditions. 
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2.8.3  Angie-of-Arrival  Fluctuations  and  ix  and  Av. 

A  key  parameter  m  determining  the  severity  of  antenna  filtering  effects  is  the 
standard  deviation  of  the  angle-of-arrival  fluctuations  Gq  of  the  electric  field  incident 
on  an  antenna.  It  is  clear  that  for  anisotropic  scattering,  when  A*  and  Ay  are  unequal, 
the  values  of  Oq  for  'ttering  about  the  x  and  y  axes  will  differ.  The  variance  of  the 
angle-of- arrival  fluctuations  about  the  x-direction  is  equal  to 
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d-Ki 

(2k, 


J  dx 
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(100) 


under  the  small-angle  scattering  approximation  that  is  required  for  the  GPSD  to  be 
valid.  A  similar  expression  holds  for  the  angle-of-arrival  variance  about  the  y- 
direction.  The  standard  deviations  of  the  angle-of-arrival  fluctuations  about  the  x  and 
y  axes  are  then  given  by 


(101a) 

(101b) 


For  small-angle  scattering  to  be  a  valid  assumption,  the  larger  of  the  two  angular 
standard  deviations  must  be  small  relative  to  one  radian.  Thus  the  minimum 
decorrelation  distance  must  be  approximately  equal  to  or  greater  than  the  RF 
wavelength  X.  The  small  angle  approximation  has  been  used  throughout  the  derivation 
of  the  GPSD,  starting  with  the  parabolic  wave  equation.  The  resulting  expression  for 
the  GPSD,  however,  does  not  exhibit  singular  behavior  w'hen  the  angle-of-arrival 
fluctuations  become  large  and  thus  the  small  angle  approximation  is  quite  robust. 

2.8.4  Isotropic  Examples. 


When  the  penetration  angle  is  zero  the  two  decorrelation  distances  are  equal  (Ax 
=  Ay  =  A())  and  the  scattering  is  isotropic  about  the  line-of-sight.  The  one-dimensional 
angle-delay  generalized  power  spectral  density  is  then  given  by 
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where  Wittwer  s  F  function  is  defined  as 
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DO 

Fiz)  =  exp  (-z;)  J  exp  (-t4-2t2z)  dt 
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Wittwer  |1980]  has  developed  a  polynomial  expansion  for  this  integral  that  is  accurate 
to  within  one  percent.  This  function  may  also  be  written  in  terms  of  K 1/4  and  1 1/4 
Bessel  functions. 

A  three-dimensional  plot  of  the  isotropic  one-dimension  angle-delay  GPSD  is 
shown  in  Figure  4.  This  plot  shows  the  mean  received  power  as  a  function  of 
normalized  angle  K£0  and  normalized  delay  coCoht.  The  vertical  axis  is  linear  with 
arbitrary  units.  The  value  of  a  is  4  for  this  figure  so  consequently  the  quantity  cocoh  is 
essentially  equal  to  27cf0. 

It  can  be  seen  that  the  power  arriving  at  large  angles  is  also  the  power  arriving 
at  long  delays.  The  power  arriving  at  long  delays  thus  has  higher  spatial  frequency 
components  than  power  arriving  at  short  delays.  When  there  is  strong  space-time 

:orrelation  (i.e.  when  \  C*t  +  Cyt  is  approximately  equal  to  unity)  these  higher  spatial 
requency  components  correspond  to  higher  Doppler  frequency  components.  The 
ignal  arriving  at  long  delays  then  varies  more  rapidly  in  time  than  the  signal  arriving 
t  short  delays. 


Figure  4.  Angle-delay  generalized  power  spectral  density. 
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Another  view  of  the  GPSD  can  be  obtained  by  considering  the  deiay-Doppler 
scattering  function 


oo 

r  dK, 


Sxd(x,(od)  =  ~2k  S(Kx,Ky,x,a)D) 


dK, 


(104a) 


For  both  the  frozen-in  and  turbulent  models  the  scattering  function  is  readily  evaluated 
for  isotropic  scattering  with  the  results: 
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(Frozen-in  Model) 


(104b) 


Std(x,cod)  =  Vrt  i0  exp 


COdXoI 

tocoh 

r  1  1 

1 

r  1 

Ct(t>cohX 

L  4  J 

2  exP 

2a2  *  wcohX 

erfc 

_V2a 

V2  J 

(Turbulent  Model) 


(104c) 


where  erfc  (•)  is  the  complementary  error  function.  It  is  clear  from  the  form  of 
Equation  104c  that  the  turbulent  model  scattering  function  is  separable  into  a  function 
of  Doppler  frequency  times  a  function  of  delay.  This  is  not  the  case  for  the  frozen-in 
model. 

For  the  general  model  the  scattering  function  can  be  written  in  the  form 
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(General  Model) 

assuming  that  Cyt  is  zero.  The  remaining  integral  in  Equation  104d  is  easily  evaluated 
numerically. 

A  comparison  of  the  scattering  functions  for  the  frozen-in  and  turbulent  models 
is  shown  in  Figure  5.  The  frozen-in  scattering  function  is  just  a  reproduction  of 
Figure  4  with  normalized  angle  Ki0  replaced  with  normalized  Doppler  frequency 
t0(od-  This  is  a  consequence  of  the  delta-function  relationship  between  angle  and 
Doppler  frequency  for  the  frozen-in  model.  For  this  model  the  signal  at  long  delays 
has  correspondingly  large  Doppler  shifts,  and  a  wing-like  structure  is  seen  in  the 
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scattering  function.  The  turbulent  model  scattering  function  does  not  exhibit  these 
Doppler  wings  because  the  Doppler  spectrum  is  the  same  at  all  delays. 

Both  functions  have  exactly  the  same  power  density  at  each  delay.  The 
difference  m  appearance  of  the  figures  is  due  to  the  fact  that  the  turbulent  model  signal 
at  long  delay  is  more  spread  out  in  Doppler  frequency  and  is  therefore  less  obvious. 

A  progression  of  scattering  functions  for  the  general  model  is  shown  in  Figure 
6.  The  space-time  correlation  coefficient  Cxt  varies  from  0.99  for  the  scattering 
function  in  the  upper  left  to  0.7  for  the  scattering  function  in  the  lower  right.  The 
scattering  function  for  Cxt  equal  to  0.99  is  essentially  identical  to  that  for  the  frozen-in 
model,  and  the  scattering  function  for  Cxt  equal  to  0.7  is  essentially  identical  to  that  for 
the  turbulent  model.  For  intermediate  values  of  CXt,  the  scattering  functions  still 
exhibit  Doppler  wings  but  the  wings  have  broader  Doppler  spectra  as  Cxt  decreases. 

2.8.5  Diffraction  Limited  Form  of  the  GPSD. 


Under  the  strong-scattering  conditions  in  the  ionosphere  that  cause  signal 
scintillation  at  radio  frequencies,  diffractive  effects  dominate  dispersive  effects.  In  this 
case  the  value  of  a  in  the  GPSD  is  large.  In  the  limit  that  a  approaches  infinity,  the 
angle-delay  part  of  the  GPSD  becomes 
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The  range  of  delay  in  this  equation  is  from  0  to  +°°  due  to  the  fact  that  the  second  term 
in  the  delta  function  is  positive. 

This  geometric  optics  limit  results  in  a  delta-function  relationship  between  angle 
and  delay.  The  delta  function  in  this  equation  can  be  used  to  associate  the  signal  with  a 
particular  angle-of-arrival  to  a  specific  delay.  This  association  is  used  to  develop 
efficient  channel  simulation  techniques. 

2.8.6  Orthogonalized  Form  of  the  GPSD  Used  in  Channel  Modelling. 

Once  the  diffraction  limit  of  the  GPSD  has  been  taken,  the  delta-function 
relationship  in  Equation  104  can  be  used  to  relate  angles-of-arrival  to  delay.  The 
GPSD  used  in  channel  modelling  can  then  be  integrated  over  delay  producing  an 
angular-Doppler  form,  denoted  Skd-  After  some  rearrangement  of  terms  this  form  of 
the  GPSD  is 
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This  fonn  of  the  GPSD  shows  that,  aside  from  the  leading  exponential  factor,  the  effect 
of  Doppler  frequency  is  to  shift  the  GPSD  in  Kx-Ky  space.  Thus  it  is  possible  to 
evaluate  Equation  106  for  zero  Doppler  frequency  and  then  to  shift  the  zero  point  of 
Kx-Ky  coordinate  system  to  obtain  the  GPSD  for  non-zero  Doppler  frequencies. 

Equation  106  does  present  a  problem  due  to  its  Kx-Ky  cross  term.  This  term 
makes  the  evaluation  of  the  signal  power  in  each  Kx-Ky  grid  cell  computationally 
difficult  and  therefore  time  consuming.  Of  course  a  simple  rotation  can  be  used  to 
eliminate  this  term.  Thus  consider  a  new  coordinate  system,  Kp-Kq,  where 

Kp  =  Kx  cosO  +  Ky  sinO  (107) 

Kq  =  -Kx  sinO  +  Ky  cosO  . 

The  rotation  angle  0  between  the  Kx-Ky  and  Kp-Kq  coordinate  system  is  chosen  to 
eliminate  the  Kp-Kq  cross  term.  This  choice  results  in  the  following  expression  for  the 
tangent  of  the  rotation  angle: 


tan  (20)  =  72 
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After  a  little  algebra  it  can  be  shown  that  the  angular-Doppler  GPSD  in  the 
Kp-Kq  coordinate  system  is 
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when  the  following  definitions  are  used: 
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Two  key  features  of  the  channel  modelling  technique  developed  for  the  general 
model  are  the  evaluation  of  the  signal  power  in  this  rotated  coordinate  system  using 
simple  error  functions,  and  the  use  of  the  Doppler  shifting  property  of  the  GPSD. 


2.9  IMPULSE  RESPONSE  FUNCTION  AND  ANTENNA  EFFECTS. 


The  channel  impulse  response  function  of  the  signal  incident  on  the  plane  of  the 
receiver  and  the  impulse  response  function  of  the  signal  at  the  output  of  an  aperture 
antenna  will  be  discussed  in  this  subsection. 

2.9.1  Channel  Impulse  Response  function. 

Consider  a  solution  U(p,o),t)  to  the  parabolic  wave  equation  in  the  plane  of  the 
receiver.  The  dependence  of  U  on  propagation  distance  z  has  been  suppressed  because 
subsequently  we  will  always  consider  U  in  the  plane  of  the  receiving  aperture.  The 
parabolic  wave  equation  solution  represents  the  random  effects  due  to  the  fluctuating 
ionosphere  on  the  incident  electric  field  at  position  p  and  time  t  from  a  transmitted 
monochromatic  wave  with  angular  frequency  (0.  The  channel  impulse  response 
function  of  the  signal  in  the  receiver  plane  is  [Knepp  and  Wittwer  1984] 


h(p,T,t) 


oo 

J  U(p,o)+(o0,t)  exp  [i<0T(co)>  +  icox] 

-OO 


dco 

2k 


(112) 


where  (0t((o))  is  the  dispersive  contribution  to  the  impulse  response  function  resulting 
from  the  total  electron  content  (TEC).  The  impulse  response  function  h(p,x,t)  is 
defined  to  be  the  received  signal  at  time  t  +  x  and  position  p  from  a  transmitted 
impulse  at  time  t  minus  the  nominal  propagation  time  (the  nominal  propagation  time  is 
usually  set  equal  to  zero). 

The  term  exp  [i(0T(O)))]  in  Equation  112  is  the  transfer  function  of  a  smooth 
ionized  plasma  and  is  equal  to  the  exponential  term  in  Equation  17.  Thus  the  smooth 
plasma  transfer  function  is 
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Because  the  smooth  plasma  or  dispersive  effects  represented  by  exp|i(0T(to))]  and  the 
fluctuating  plasma  effects  represented  by  U(p,zr,co,t)  appear  as  the  product,  it  is 
possible  to  separately  consider  these  two  effects.  The  dispersive  effects  will  be 
considered  in  the  next  subsection. 


If  the  transmitted  signal  is  a  modulated  waveform  mt(t),  then  the  signal  complex 
voltage  incident  on  the  plane  of  the  receiver  is  the  convolution  of  the  transmitted 
modulation  and  the  channel  impulse  response  function: 


oo 

u(p,t)  =  J  mt(t-T)  h(p,T,t-T)  dr  . 

-  oo 


(114) 


It  is  clear  from  the  form  of  this  equation  that  received  signal  at  time  t  and  delay  t 
corresponds  to  transmitted  modulation  at  time  t  -  T.  The  third  argument  of  h  in  the 
equation  is  t  -  I  because  modulation  transmitted  at  earlier  times  transits  the  channel  at 
earlier  times.  It  is  generally  assumed  that  the  channel  varies  slowly  in  time  relative  to 
the  delay  spread  of  the  signal,  and  h(p,T,t-T)  is  set  equal  to  h(p,T,t)  in  Equation  114. 
This  is  strictly  true  only  if  the  product  f0T0  is  large.  However,  the  product  f0T0  will 
generally  be  large  under  conditions  where  RF  signals  can  propagate  through  the 
ionosphere. 

If  the  delay  spread  of  the  impulse  response  function  is  larger  than  the  modulation 
period  of  mt(t),  then  the  convolution  in  Equation  114  will  encompass  multiple 
modulation  periods.  As  a  result  the  received  signal  will  contain  information  from 
multiple  modulation  periods.  In  communications  systems  this  phenomenon  is  referred 
to  as  intersymbol  interference. 

2.9.2  Dispersive  Effects. 

When  the  dispersive  term  (0j(co))  is  expanded  in  a  Taylor  series  about  the 
carrier  radian  frequency  the  result  is 

2 

(0T(C»)  =  (0T(O3o))  -  (C0-0)o)  (0T(O)O)> +^Y^-<0T(O)O)>  +  ••  (115) 

where  the  first  three  coefficients  in  the  expansion  are 
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These  equations  may  be  expanded  using  the  assumption  that  the  carrier  frequency  is 
much  larger  than  the  plasma  frequency  with  the  results: 


(0T((%)>  =  -  +  ^rcNx 

/o',  R  ^cNt 

<0T(wo)>  =  -+~2sr 
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where  the  free  space  range  R  and  total  electron  content  Nj  are  given  by 
R  =  zt  +  zr 
and 
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The  first  terms  in  Equations  117a  and  117b  are  simply  the  free  space  phase  shift  and 
propagation  time  which  are  proportional  to  the  line-of-sight  distance  R.  The  terms 
proportional  to  Nj  in  Equations  117  represent  the  phase  shift,  group  delay  and 
dispersion  due  to  the  mean  ionization. 
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The  Doppler  shift  fo  of  the  incident  signal  due  to  range  and  TEC  dynamics  is 


J_  d(0i(coo))  __  1  dR  Xrc  dNj 

tD  ~  2 k  dt  ~  '  k  dt  +  2k  dt 


(120) 


Note  that  increasing  TEC  (positive  dNx/dt)  increases  both  the  propagation  time  and  the 
Doppler  shift  whereas  increasing  R  (positive  dR/dt)  increases  the  propagation  time  but 
decreases  the  Doppler  shift.  » 

2.9.3  Antenna  Aperture  Effects. 


The  voltage  at  the  output  of  an  aperture  antenna  is  the  spatial  convolution  of  the 
incident  voltage  and  the  aperture  weighting  function,  A(p).  The  received  voltage  for 
an  antenna  located  at  p0  and  pointing  in  the  K0  direction  is  given  by 


OO 

UA(po,(o,t)  =  JU(p,co,t)  A(po-p)  exp  (iKo  (po-p)]  d2p  (121) 

-CO 


where  the  subscript  A  denotes  the  voltage  at  the  output  of  the  antenna.  The  z 
dependence  of  Ua  has  been  suppressed  because  it  is  understood  that  this  voltage  is  at 
the  plane  of  the  receiver.  It  is  assumed  that  the  aperture  weighting  function  is 
independent  of  frequency.  This  is  generally  true  for  a  range  of  frequencies  about  the 
carrier  frequency  that  is  larger  than  the  signal  bandwidth. 

In  order  to  relate  the  GPSD  of  Ua  to  the  GPSD  of  the  incident  signal,  the  two- 
position,  two-frequency,  two-time  mutual  coherence  function  of  Ua  is  required.  The 
mutual  coherence  function  of  the  signal  out  of  the  antenna  is 

TA(p,0),t)  =  <UA(Pi,0)1,ti)UA(p2.W2.t2))  (122) 
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x  A(p'-p1)A*(p"-p2)  exp  [iKo-CPi  -  p'  -  p2  +  P")l  • 

For  statistically  stationary  processes  the  expectation  of  UaUa  must  be  a  function  only 
of  the  differences  p  =  pj  -  p2.  co  =  C0|  -  co2,  and  t  =  t]  - 12,  and  the  expectation  of  UU* 
in  the  integrand  must  be  a  function  only  of  the  differences  p'  -  p”,  to,  and  t. 

The  aperture  weighting  function  can  be  written  in  terms  of  the  voltage  beam 
pattern  g(Kj_)  using  the  Fourier  transform  relationship 


r  d2Ki 

A(p)  =  J  g(Ki)  exp  (iKi-p)  . 


023) 


38 


Upon  substituting  this  equation  for  both  aperture  weighting  functions  in  the  expression 
for  the  mutual  coherence  function,  changing  variables  from  p'  to  r  =  p’  -  p",  and 
changing  the  order  of  integration.  Equation  122  becomes 
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rA(p,co,t)  =  exp  (iK0  p)  J  d2r  r(r,co,t)  exp  (-iK0-r) 
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The  last  integral  in  this  expression  is  equal  to  (2k)2  5(K"-K’),  and  the  K"  integral  may 
be  performed  directly.  Another  change  in  the  order  of  integration  yields 

oo  2 

TA(p,co,t)  =  exp  (iK0  p)  J  G(K')  exp  (iK'-p)  (125) 

oo 

x  J  d2r  T(r  ,co,t)  exp  [-i(K'  +  K0)r]  . 

-oo 

The  quantity 

G(K)  =  g(K)g*(K)  (126) 

is  the  power  beam  pattern  of  the  antenna. 

The  mutual  coherence  function  r(r,o),t)  of  the  signal  incident  on  the  plane  of  the 
antenna  that  appears  in  the  second  integral  of  Equation  125  is  the  product  of  the  free 
space  term  T0  (Eqn.  43)  and  the  stochastic  term  Tj.  The  free  space  term  may  be  pulled 
out  of  the  second  integral  if  it  is  assumed  not  to  vary  over  the  face  of  the  antenna.  This 
is  equivalent  to  assuming  that  any  deviations  from  a  plane  wave  in  the  incident  signal 
are  due  to  scattering  effects  in  the  ionosphere  and  are  not  due  to  geometrical  effects. 
After  the  free  space  term  is  pulled  out  of  the  integral,  Ta  may  be  assumed  to  represent 
only  the  stochastic  fluctuations  of  the  received  signal. 

Now  the  GPSD  of  the  signal  out  of  the  antenna  can  be  computed  by  taking  the 
appropriate  Fourier  transforms  (see  Equation  84)  from  p,  o),  and  t  to  Kj_,  x,  and  cod 
respectively.  After  performing  the  co  to  x  and  t  to  cod  transforms  and  rearranging 
terms,  the  GPSD  at  the  antenna  output  is 
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The  last  integral  in  this  equation  is  just  (2 kY  8(K'-Kj_+K0).  Using  this  delta  function 
to  perform  the  K'  integral  results  in: 


Sa(Ki,t,coD)  =  G(Ki-K0)  j  d2r  f(r,x,(OD)  exp  (-iK±r)  (128) 

-co 

where  the  remaining  integral  in  this  expression  is  just  S(Ki,t,cod).  The  GPSD  of  the 
signal  out  of  the  antenna  is  then 

Sa(Kx,t,cod)  =  G(Ki-K0)  S(Ki,t,cod)  ■  129) 

As  expected,  the  effect  of  an  antenna  is  to  modify  the  incident  power  as  a 
function  of  angle.  This  result  will  be  used  throughout  the  rest  of  this  report. 

The  effect  of  antenna  filtering  is  illustrated  in  Figure  7  which  shows  four  plots 
of  the  GPSD  at  the  outputs  of  four  uniformly-weighted  circular  antennas  (with 
diameter  D)  for  isotropic  scattering  (ix  =  5.y  =  &o).  The  antennas  in  this  example  are 
all  pointed  along  the  line-of-sight.  The  plot  in  the  upper  left  is  the  same  as  in  Figure  4 
and  is  for  a  point  antenna  (D  «  io),  f°r  which  there  is  no  antenna  filtering.  The 
other  three  plots  in  the  figure  are  for  cases  where  the  ratio  of  the  antenna  diameter  to 
the  decorrelation  distance  is  equal  to  1,2,  and  4.  As  the  antenna  size  increases  for  a 
given  value  of  K0  or,  equivalently,  as  the  decorrelation  distance  decreases  for  a  given 
antenna  diameter,  more  of  the  incident  signal  with  large  angles-of-arrival  is  filtered 
out  of  the  received  signal.  The  effect  of  this  on  the  mean  received  power,  frequency 
selective  bandwidth,  and  decorrelation  time  of  the  signal  out  of  the  antenna  will  be 
discussed  next. 
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SECTION  3 

ANTENNA  FILTERING  EFFECTS 

An  antenna  beam  pattern  acts  as  an  angulai  filter  of  the  received  signal  power. 
Because  of  this,  the  mean  power,  decorrelation  time,  and  frequency  selective 
bandwidth  of  the  signal  at  the  output  of  an  antenna  are  different  than  those  of  the 
incident  signal.  The  reduction  in  mean  power  is  a  direct  consequence  of  the  attenuation 
of  the  signal  arriving  outside  of  the  main  beam  of  the  antenna.  If  the  antenna  is 
pointed  along  the  line-of-sight,  the  signal  that  is  attenuated  by  the  antenna  is  that  with 
large  angles-of-arrival  and  hence  large  times-of- arrival.  Thus  the  effect  of  ai.  antenna 
pointed  along  the  line-of-sight  is  to  reduce  the  delay  spread  of  the  signal.  When  the 
same  antenna  is  pointed  away  from  the  line-of-sight,  however,  its  effect  on  the  delay 
spread  of  the  output  signal  can  be  quite  different,  as  will  be  shown  in  this  section. 

Depending  on  the  values  of  the  correlation  coefficients  Cxt  and  Cyt,  the  antenna 
may  also  filter  the  decorrelation  time.  In  the  fully  turbulent  case,  the  signal  at  all 
angles-of-arrival  has  the  same  effective  decorrelation  time,  and  hence  an  antenna  has 
no  effect  on  the  decorrelation  time  of  the  output  signal.  However  in  the  frozen-in  case 
the  effective  decorrelation  time  of  the  signal  varies  with  angle-of-arrival,  and  an 
antenna  can  strongly  filter  the  decorrelation  time.  This  filtering  results  in  a  signal  that 
is  more  slowly  varying  at  the  antenna  output  than  it  is  on  the  face  of  the  antenna. 

The  effects  of  aperture  antennas  with  arbitrary  beamwidths  and  pointing  angles 
will  be  considered  in  this  section.  The  antenna  beam  patterns  for  uniformly-weighted 
circular  or  rectangular  apertures  and  for  Gaussian  apertures  are  described  in  Section 
3.1.  Then  the  filtering  equations  for  mean  power,  spatial  and  temporal  decorrelation 
properties,  and  frequency  selective  bandwidth  of  the  signal  out  of  a  Gaussian  antenna 
with  arbitrary  beamwidths  and  pointing  angles  are  given  in  Section  3.2.  Many  of  these 
results  have  been  derived  previously  by  Frasier  [1988]  and  are  just  reproduced  here. 

An  example  is  given  in  Section  3.3.  Here  the  antenna  is  assumed  to  be  circular 
with  uniform  weighting,  and  the  scattering  is  assumed  to  be  isotropic  (ix  =  £y  =  £0). 
The  filtering  equations  for  this  case  are  plotted  versus  the  ratio  of  the  antenna  diameter 
to  the  decorrelation  distance  for  a  range  of  pointing  angles.  This  example  illustrates 
most  of  the  important  effects  of  antenna  filtering. 

3.1  ANTENNA  DESCRIPTIONS. 

The  coordinate  systems  for  the  propagation  and  the  antenna  are  shown  in  Figure 
8.  The  z  axis  is  along  the  line-of-sight,  and  it  is  assumed  that  the  face  of  the  antenna  is 
in  the  x-y  plane.  Note  that  for  a  receiving  antenna  we  have  reversed  the  direction  of 
the  z  axis  in  the  figure  from  pointing  in  the  direction  of  propagation  to  pointing  back 
toward  the  transmitter.  This  is  done  because  it  is  usual  for  antenna  coordinate  systems 
to  point  out  of  rather  than  into  an  aperture.  The  beam  of  the  antenna  in  the  figure  is 
pointed  away  from  the  line-of-sight  in  the  direction  K0  defined  by  a  pointing  direction 
elevation  angle  ©0  (measured  from  the  z  axis)  and  an  azimuth  angle  O0  (measured 
from  the  u  axis).  The  rotation  angle  T  is  the  angle  between  the  scattering  x  axis  and 
the  antenna  u  axis. 
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Figure  8.  Propagation  and  antenna  coordinate  systems. 

The  antenna  (u-v)  coordinate  system,  where  v  is  in  the  x-y  plane  and  is 
orthogonal  to  u,  is  chosen  for  convenience  in  describing  the  antenna  beam  pattern.  For 
example,  if  the  antenna  is  rectangular,  the  (u-v)  axes  should  be  aligned  with  the  sides  of 
the  aperture.  It  is  assumed  that  the  antenna  beam  pattern  is  separable  in  (u-v) 
coordinates.  Thus  the  two-dimensional  beam  pattern  G(0U,9V)  is  approximated  by  the 
product  G(0u)G(0v)  where  9U  and  0V  are  angles  about  the  u  and  v  axes  respectively.  In 
many  cases  the  u  axis  will  be  parallel  to  the  local  earth  tangent  plane. 


3.1.1  Gaussian  Beam  Profiles. 


The  antenna  beam  pattern  is  assumed  to  be  separable  in  the  u-v  coordinate 
system.  The  Gaussian  antenna  beam  pattern  is  then 

G(Ku,Kv)  =  exp  au(Ku-Kou)  -  cXy(Kv-Kov)  J  (130) 


where  Kou  and  Kov  are  components  of  the  pointing  vector.  The  peak  gain  G(0,0)  in 
this  equation  has  been  set  to  unity  because  this  value  is  usually  included  in  the 
calculation  of  the  mean  received  power.  The  components  of  the  pointing  vector  in  the 
u-v  plane,  K()u  and  Kov,  are: 
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For  either  the  u  or  the  v  direction,  the  antenna  beam  pattern  for  pointing  along 
the  line-of-sight  can  also  be  written  as 


' 

(2Qf 

exp 

-  In2 

l0oj , 

(132) 


where  90  is  the  3-dB  beamwidth  (i.e.,  full  width  at  half  maximum).  Equating  the  two 
patterns  gives 


au 


1  n 2  X2 

rc2e02u 


(133) 


where  90u  is  the  3-dB  beamwidth  in  the  u-direction.  A  similar  expression  holds  for  a2 
which  relates  this  parameter  to  the  3-dB  beamwidth  in  the  v-direction,  90v. 


3.1.2  Uniformly  Weighted  Circular  Apertures. 

For  a  uniformly-weighted  circular  aperture,  the  aperture  weighting  function  is 


A(p) 


'I 
<  a 

,0 


if  Ipl  <  y 
otherwise 


(134) 


where  D  is  the  diameter  of  the  circular  aperture.  The  value  of  a  is  chosen  so  that  the 
peak  antenna  gain  is  unity.  The  voltage  antenna  gain  pattern  g(KjJ  is  related  to  the 
aperture  distribution  function  by  the  Fourier  transform 
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g(Kj.)  =  J  exp  (-iKi  p)  A(p)  d2p 


(135) 


With  p  confined  to  the  plane  of  the  aperture,  the  dot  product  K±  p  is 
Kj_  p  =  kr  sin  8  sin  <p 


(136) 


where  0  is  the  angle  from  the  normal  to  the  aperture,  <p  is  the  azimuth  angle,  k  =  2k/X, 
and  r  =  Ipl.  The  Fourier  transform  then  becomes 


g(9)  =  —  J  d  r  r  J  dcp  exp  (-ikr  sin  8  sin  cp)  . 


(137) 


0  -k 


Performing  the  indicated  integrals  results  in  the  well  known  form  for  the  power  beam 
pattern 


G(8)  =  g(0)g*(0)  = 


4J2  [7t(DA)sin81 
[TC(DA)sin8]2 


(138) 


where  Jt  is  the  Bessel  function,  and  where  G(0)  =  1  when  the  value  of  a  is  chosen  to  be 
the  area  of  the  aperture,  ttD2/4. 

The  3-dB  (full  width  at  half  maximum)  beamwidth  [i.e.  G(0o/2)  =  1/2]  is  given 
by  solving  the  equation 


4Jf(x)  _  l 
x2  ~  2 


0  39) 


with  the  result  x  =  1.616340.  Assuming  that  the  beamwidth  is  small  so  sin(0o/2)  « 
Qq/2,  the  beamwidth  in  terms  of  the  diameter  D  is 


0O  =  1.02899  ^  radians 


(140) 


If  the  uniformly-weighted  circular  antenna  beam  pattern  is  approximated  by  a  Gaussian 
pattern  with  the  same  3-dB  beamwidth,  then  the  a  coefficients  that  appear  in  Equation 
129  are 


„2  _  „2  _  _l!n2]Di_ 
u  v  (1.02899TC)2 


(141) 
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3.1.3  Uniformly  Weighted  Rectangular  Apertures. 

For  uniformly-weighted  rectangular  apertures,  the  aperture  weighting  function 
is 


A(u,v)  =  < 


Du  ,, 

<  ~y  and 


Dy 
-  2 


(142) 


10  otherwise 

where  Du  and  Dv  are  the  lengths  of  the  aperture  in  the  u  and  v  directions  respectively. 
In  this  case,  the  Fourier  transform  indicated  in  Equation  1 35  gives  the  result 

G(0u,ev)  =  G(0U)  G(6V)  (143) 


where  G(0U)  has  the  familiar  sin2(x)/x2  form: 

sin2  [7t(Du/V)sin0E] 

G(e“>  ■  ' 


(144) 


A  similar  expression  holds  for  G(0V).  The  normalization  a  of  the  aperture  weighting 
function  is  just  the  area  DUDV  of  the  rectangular  antenna.  The  3-dB  beamwidth  is  given 
by  solving  the  equation 


sin2(x)  _ 
x2  ~  2 


(145) 


with  the  result  x  =  1.391557.  The  u  and  v  beamwidths  in  terms  of  the  antenna  sizes  Du 
and  Dv  are  then  given  by 


0OU  =  0.885893 


Du 


0, 


0v 


-  0.885893 


X_ 

Dv 


(radians)  (146) 


assuming  that  sin(0o/2)  =  0o/2.  If  this  beam  pattern  is  approximated  by  a  Gaussian 
pattern  with  the  same  3-dB  beamwidth,  the  a  coefficients  that  appear  in  Equation  1 30 
are 

2  (l"2)Du  2  C"2)P? 

““  (0.885893IC)2  “v  (0.885893JI)2  ' 

3.2  FILTERING  EQUATIONS. 

The  filtering  equations  relate  the  statistics  of  the  signal  at  the  outputs  of  one  or 
more  antennas  to  the  statistics  of  the  signal  incident  on  the  antennas.  The  statistics  that 
are  considered  in  this  section  are  the  mean  power,  frequency  selective  bandwidth, 
decorrelation  distances,  and  decorrelation  time  of  the  signal  out  of  an  antenna  and  the 
cross  correlation  of  the  signals  out  of  two  separate  antennas.  Gaussian  beam  patterns 
will  be  assumed  for  mathematical  convenience.  An  elegant  derivation  of  most  of  the 
equations  given  in  this  section  may  be  found  in  Frasier  [1988]. 
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The  severity  of  the  filtering  effects  is  determined  by  the  relative  size  of  the 
standard  deviation  of  the  angle-of-arrival  fluctuations  ae  and  the  antenna  beamwidths. 
When  oq  is  small  compared  to  the  antenna  beamwidths,  the  signal  arrives  essentially  at 
the  peak  of  the  beam  pattern,  if  the  pointing  error  is  small,  and  the  filtering  effects  are 
small.  However,  if  <79  is  large  compared  to  the  beamwidths,  much  of  the  signal  arrives 
at  angles  outside  of  the  main  lobe  of  the  antenna  beam  pattern  and  the  filtering  effects 
are  large.  Equivalently,  large  values  of  the  ratio  ae/0o  correspond  to  situations  where 
the  decorrelation  distance  of  the  incident  signal  is  small  compared  to  the  antenna  size 
and  the  incident  electric  field  as  seen  by  the  aperture  is  no  longer  a  plane  wave.  In  this 
situation,  the  induced  voltages  across  the  face  of  the  aperture  do  not  add  coherently 
when  summed  together  by  the  antenna  with  a  loss  in  signal  power  as  a  result. 

In  the  next  subsections,  expressions  will  be  presented  for  the  mean  power, 
frequency  selective  bandwidth,  decorrelation  time,  and  decorrelation  distances  of  the 
signal  at  the  output  of  an  antenna.  Clearly  the  numerical  value  of  these  quantities  must 
be  independent  of  the  choice  of  coordinate  systems  (e.g.,  x-y  or  u-v  coordinates)  and 
Frasier  [1988]  has  derived  expression  for  these  quantities  which  are  coordinate-system 
independent.  However  evaluation  of  the  filtering  equations  requires  a  choice  of 
coordinate  systems,  and  the  x-y  systems  will  be  used  here. 


Before  presenting  the  filtering  equations,  it  is  convenient  to  define  some 
quantities  that  are  common  to  all  these  equations: 


4otu  ,  4al  , 

Qx  -  1  +  02  cos2vF  +  2  sin2vF 

(148) 

4  u.1  ,  4otv  , 

Qy  =  1  +  2  s>n2vE  +  2  costVE 
*  y  *y 

(149) 

4(otu  -  a\) 

Qxy  =  cosH'sinH' 

(150) 

Qo  =  V QxQy  -  Q xy  • 

(151) 

The  components  of  the  pointing  vector  in  the 
and  are  obtained  from  their  form  in  the  u-v 
simple  rotation: 

x-y  coordinate  system  are  also  required 
coordinate  system  (Eqn.  131)  using  a 

K0x  =  Y  sine°  cos  (°o+vf/) 

(152a) 

2k 

Koy  =  y  sin©o  sin  (dV^)  • 

(152b) 
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3.2.1  Scattering  Loss. 

The  mean  power  of  the  signal  out  of  an  antenna  is  calculated  using  the 
expression: 


Pa  =  j 


7d2Ki 


OO  CO 


JdT  fd-^f  G(K±- K0)  S(K±fT,<BD) 


(153) 


-OO  -CO 


which  gives  the  result 


Pa  =  n  exP  '  1  ' 


Qy  \  K  Qx 


'■  Q  l 


0 'jl)  Qxy 

Qo  j  4  '  2Q2KoxUKoyay  . 


0  54) 

The  first  term  in  this  equation,  1/Qo,  is  just  the  mean  received  power  when  the  antenna 
is  pointed  along  the  line-of-sight. 

In  the  limit  that  both  decorrelation  distances  are  large  relative  to  the  size  of  the 
antenna  only  the  exponential  term  in  the  expression  for  Pa  differs  from  unity,  and  the 
expression  for  the  mean  received  power  reduces  to 


Limit 

4X  _  oo  ?A  =  exp  [.(a2K2u  +  o^j^j 
!y  -»  oo 

which  is  equal  to  the  antenna  beam  gain  in  the  direction  of  the  line-of-sight. 
The  scattering  loss  of  the  antenna  in  decibels  is 
LS(dB)  =  -10  logio  (PA)  • 


(155) 


(156) 


This  loss  therefore  includes  both  the  loss  in  the  mean  received  power  due  to 
scintillation  and  the  loss  due  to  the  fact  that  the  antenna  may  be  pointed  away  from  the 
line-of-sight. 

3.2.2  Frequency  Selective  Bandwidth. 

The  frequency  selective  bandwidth  is  defined  in  Section  2.8.2  in  terms  of  the 
time  delay  jitter  of  the  signal.  At  the  output  of  an  antenna,  these  moments  are  given  by 


7d2K! 

PA<tn>  -  J 


oo  oo 

J  dx  xn  j  “  G(Ki-  K0)  S(Ki,t,(0D)  . 


(157) 


Straightforward  evaluation  of  the  indicated  integrals  for  n  equal  to  1  and  2  is  indeed  a 
formidable  algebraic  task.  However  Frasier  [1988]  gives  surprisingly  simple¬ 
appearing  expressions  for  at  and  the  antenna-filtered  value  of  the  frequency  selective 
bandwidth,  fA,  which  are  independent  of  the  choice  of  the  coordinate  system.  In  terms 


of  the  notation  used  in  this  report,  the  expression  for  fA  is  somewhat  more  complicated 
and  is 


2 

q  2  n  2  Qo 

.  xxxy  . _ 


(158) 


*xQx  HQy  2 
~rx  +  +  2Q  xy 

.  xy  xx 


where  the  term  P  gives  the  effects  of  antenna  pointing.  This  pointing  angle  term  may 
be  written  as: 


P  -  PxK02xHx  +  PyKo^lJ  +  PxyKoxlxKoyly 


where 


Px  =  2QxV(1-Qx)  + 


Py  =  2Q2xy(l-Qy)  + 


Qy  ®-y(Qo'Qy)2  QxQxyM 


Qxix(Qo-Qx)2  QyQxyi( 


Pxy  =  Qxy  (Qx-D(l-Qy)-Qxy  + 


3.2.3  Decorrelation  Time. 


Q*a2*(Qo-Q*)  Qy  ^y(Qo'Qy) 


«  2  + 


(159) 


(160a) 


(160b) 


(160c) 


The  temporal  coherence  function  of  the  signal  out  of  an  antenna  is  given  by 


r  via*  I  r  r  uuj  I  i 

PArA(0  =  j  7Z~2  /  dT  /  -2^exp(-iC0Dt)SA(Kx,Ky,T,C0D)  .  (161) 

-oo  '  '  -oo  -oo 


In  general,  the  temporal  coherence  function  is  complex  when  the  pointing  angle  is  non¬ 
zero  because  antenna  pointing  results  in  a  mean  Doppler  shift  of  the  output  signal.  The 
antenna-filtered  decorrelation  time  ta  is  then  calculated  by  finding  the  1/e  point  of 
1  Ta  ( t) I  with  the  result 


2  n2  CxlQy  +  CytQx  -  2CxtCylQxy  2 
xt  yl+  Qo 


In  the  turbulent  limit  (i.e.,  Cxt  =  0  and  Cyt  =  0), 


(162) 


(Turbulent  Model) 


(163) 
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and  in  the  frozen-in  limit  (i.e.,  Cxt  =  1  and  Cyt  =  0) 

(Frozen-in  Model)  (164) 

It  is  interesting  to  note  that  the  antenna-filtered  decorrelation  time  is  not 
explicitly  dependent  on  pointing  angle,  but  does  depend  on  the  ratio  of  the  antenna 
beamwidth  and  the  standard  deviation  of  the  angular  scattering  through  the  Q  factors. 
Of  course  as  a  phased  array  antenna  beam  is  scanned  away  from  boresight,  the 
effective  aperture  size  will  decrease  thereby  broadening  the  beamwidth  and  implicitly 
changing  the  value  of  the  filtered  decorrelation  time. 

The  mean  Doppler  shift  (Oa  due  to  antenna  pointing  is 

[Cxt(Qo'Qy)'CytQxy]Kox  fi.x  -  [Cyt(Qo-Qx)-CxtQxy]Koyiy 

“A  =  t„Q l  ■  0  65) 

This  quantity  is  clearly  equal  to  zero  for  the  turbulent  model.  For  non-zero  values  of 
the  space-time  correlation  coefficients,  the  mean  Doppler  shift  is  proportional  to  the 
components  of  the  pointing  vector. 

3.2.4  Decorrelation  Distances. 

The  x-direction  decorrelation  distance  of  the  signal  out  of  the  antenna  is  given  by 
the  1/e  point  of  Ta(x)  where 

oo  oo  oo  oo 

,dKv  r  r  dcon  r  dKx 

PArA(x)  =  J  2k  j  dr  J  J  *2^- exp  (iKxx)  SA(Kx,Ky,T,coD)  .  (166) 

-  oo  -oo  -oo  -oo 

A  similar  expression  holds  for  TA(y).  These  two  expression  give  the  following  results 
for  the  decorrelation  distances  i  Ax  and  5.  Ay  of  the  signal  at  the  antenna  output: 


TA  = 

T°  Vq7  ‘ 


2- Ax 

Qo 

(167) 

VQy 

2  Ay 

Qo 

2y 

V  Qx  ' 

(168) 

The 

interpretation  of  these  quantities  is  that  2ax 

and  4  Ay  are  the  distances  in  the 

x-  or  y-direction  respectively  that  the  antenna  must  be  instantaneously  displaced  in 
order  for  the  normalized  cross  correlation  of  the  output  signal  to  have  a  value  of  1/e. 
As  was  the  case  for  the  antenna-filtered  decorrelation  time,  these  quantities  do  not 
explicitly  depend  on  pointing  angle. 
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3.3  ISOTROPIC  EXAMPLE. 


An  isotropic  example  is  presented  in  this  section  to  illustrate  some  of  the  effects 
of  the  antenna  on  the  parameters  discussed  above.  For  this  example,  both  the  angular 
scattering  and  the  antenna  beam  pattern  will  be  assumed  to  be  isotropic,  so 


iX  -  iy  -  i0 

and 

2  ,  ln2  X2 

=  °v  =  W 


069) 


(170) 


Without  further  loss  of  generality,  it  can  be  assumed  that  the  antenna  is  pointed  in  the 
x-direction  and  that  the  rotation  angle  is  zero.  The  components  of  the  pointing  vector 
are  then  given  by 


Kox 


2k 

X 


sin  0O 


(171a) 


K  0y  —  0  .  (171b) 

With  these  assumptions  the  Q  factors  become 

Q,.qy-Q«,+4^g  (172a) 

Qxy  =  0  .  (172b) 

At  this  point  it  is  convenient  to  write  the  antenna  beamwidth  in  terms  of  the  antenna 
size  in  order  to  eliminate  the  explicit  dependence  of  the  Q  factors  on  carrier 
wavelength: 


Go  -  ao 


X 

D 


(173) 


where,  for  uniformly-weighted  antennas, 

f  1.02899  Circular  antenna 

I  0.885893  Square  antenna 

The  Q  factor  can  then  be  written  as 


D 

1  +  0.265  -2 


D 

1  +  0.358  -2 


Circular  antenna 


Square  antenna 


(174) 


(175) 
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The  mean  received  power  for  isotropic  scattering  and  an  isotropic  antenna  is 


Pa 


1 

Q  exP 


aotr2(Q-l)  ©o 

Q  06 


(176) 


when  the  pointing  angle  is  assumed  to  be  small  compared  to  one  radian  so  sin@0  is 
approximately  equal  to  ©0.  For  a  given  value  of  &(yD  <\c.,  for  a  fixed  value  of  Q)  the 
effect  of  antenna  pointing,  as  expected,  is  to  mono^mcally  decrease  the  mean  received 
power  as  the  antenna  beam  is  pointed  farther  ..vay  from  the  line-of-sight.  In  the  limit 
that  the  decorrelation  distance  is  much  |,'.,ger  than  the  antenna  diameter,  the  Q  factor 
approaches  unity  and  the  mean  ’•.ceived  power  is  approximately  equal  to  1/Q 
independent  of  pointing  angle  in  other  words,  when  the  angular  scattering  standard 
deviation  is  large  compared  to  the  antenna  beamwidth  and  the  pointing  angle,  the  mean 
power  is  insensitive  to  me  location  of  the  beam  peak  relative  to  the  line-of-sight. 

The  scattering  loss  L$  (L$  =  1/Pa)  for  a  uniformly-weighted  circular  antenna  is 
plotted  in  Figure  9  versus  the  ratio  £(/ D  when  the  pointing  angle  of  the  antenna  beam 
is  0,  0o/2  and  0O.  It  is  assumed  for  this  example  that  the  beamwidth  remains  constant  as 
the  beam  in  pointed  away  from  the  line-of-sight.  When  the  decorrelation  distance  is 
large  compared  to  the  antenna  diameter,  the  curves  approach  the  line-of-sight  loss  of 
beam  (0  dB  for  a  pointing  angle  equal  to  zero,  3  dB  for  a  pointing  angle  equal  to  one- 
half  beamwidth,  and  12  dB  for  a  pointing  angle  equal  to  a  full  beamwidth).  When  the 
decorrelation  distance  is  small  compared  to  the  antenna  diameter,  the  angle-of-arrival 
spread  of  the  incident  signal  is  large  compared  to  the  beamwidth,  so,  as  mentioned 
above,  the  scattering  loss  is  insensitive  to  pointing  angle  as  long  as  the  pointing  angle 
remains  within  a  few  beamwidths  of  the  line-of-sight.  Between  these  two  limits,  the 
scattering  loss  of  a  beam  pointed  away  from  the  line-of-sight  may  actually  decrease 
with  decreasing  Hq/D  as  the  signal  is  scattered  away  from  the  line-of-sight  and  into  the 
main  lobe  of  the  antenna. 

The  equation  for  the  ratio  of  the  filtered  to  unfiltered  frequency  selective 
bandwidth,  fA/fo.  for  this  isotropic  example  is 


fA 

fo 


Q  +  2a207t2(Q-l) 


2©0  lo 
0o  D2 


i_ 

2 


(177) 


It  is  interesting  to  note  that  when  the  pointing  angle  is  zero  the  ratio  fA/fo  is  equal  to  Q 
which  is  also  equal  to  the  scattering  loss  in  this  case.  Figure  10  shows  plots  of  the  ratio 
fA/fo  as  a  function  of  the  ratio  io/D  for  a  uniformly-weighted  circular  antenna  and  for 
three  pointing  angles. 
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Figure  9.  Scattering  loss  for  a  uniformly-weighted  circular  antenna  and 
isotropic  scattering. 


lo/D 

Figure  10.  Filtered  frequency  selective  bandwidth  for  a  uniformly- 
weighted  circular  antenna  and  isotropic  scattering. 
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For  an  antenna  pointed  along  the  line-of-sight,  the  potentially  beneficial  effect  of 
antenna  filtering  is  that  the  signal  which  is  filtered  out  is  that  which  arrives  at 
relatively  large  delays.  It  is  this  delayed  signal  that  causes  most  of  the  intersymbol 
interference  in  the  detection  and  demodulation  of  wide  bandwidth  signals.  Stated 
another  way,  because  the  antenna  filters  out  much  of  the  delayed  sir  at,  and  since  the 
frequency  selective  bandwidth  is  an  inverse  measure  of  the  signal  delay  spread,  the 
filtering  increases  the  frequency  selective  bandwidth  of  the  output  signal  relative  to  that 
at  the  antenna  input. 

Clearly  the  situation  is  different  if  the  antenna  is  pointed  away  from  the  line-of- 
sight.  The  ratio  fA/fo  is  less  than  unity  for  non-zero  pointing  angles  and  values  of  Hq/D 
between  about  0.1  and  1.0.  This  implies  that  the  standard  deviation  of  the  delay  jitter 
is  increased  by  the  antenna.  This  does  not  imply  however  that  the  antenna  is  somehow 
creating  more  signal  power  at  long  delays  than  is  incident  on  the  antenna,  as  measured 
by  f0.  Rather,  the  antenna  pointed  away  from  the  line-of-sight  is  weighting  the  power 
at  long  delays  more  than  that  at  short  delays  with  the  peak  weighting  occurring  at  the 
delay  that  corresponds  to  the  pointing  angle.  The  result  is  an  increased  value  of  the 
time  delay  jitter  standard  deviation.  Thus  intersymbol  interference  effects  in  the 
output  of  an  antenna  beam  may  increase  as  the  beam  is  pointed  away  from  the  line-of- 
sight. 


Antenna  filtering  also  affects  the  decorrelation  time  of  the  received  signal.  For 
isotropic  scattering  and  an  isotropic  antenna,  the  equation  for  the  ratio  Ta/Tq  reduces  to 


T°  V Q  +  (C2,i+C£,)(l-Q)  ' 

This  equation  is  plotted  in  Figure  11  versus  the  ratio  Hq/D  for  various  values  of 

C  =  VcL  +  C2yt  .  (179) 

As  discussed  above,  the  antenna-filtered  value  of  decorrelation  time  is  not 
explicitly  dependent  on  the  pointing  angle.  It  does  depend  strongly  on  the  model  for 
the  temporal  fluctuations.  For  the  turbulent  model  (C^i  +  Cyt  =  0),  the  Doppler 
spectrum  is  independent  of  angle,  so  the  signal  arriving  at  all  angles  has  the  same 
decorrelation  time.  Hence  the  filtered  value  of  decorrelation  time  is  equal  to  that  of 
the  incident  signal  (i.e.,  taAo  =  1)  for  this  model.  For  larger  values  of  Cxt  +  Cyt  there 
is  coupling  between  the  angular  spectrum  and  the  Doppler  spectrum  of  the  incident 
signal.  The  effect  of  an  antenna  is  to  narrow  the  angular  spectrum  of  the  received 
signal  as  shown  in  Figure  7.  Thus  when  +  Cyt  is  greater  than  zero  the  antenna  also 
narrows  the  Doppler  spectrum  of  the  received  signal.  Because  the  decorrelation  time 
of  the  signal  is  an  inverse  measure  of  the  width  of  the  Doppler  spectrum,  the 
decorrelation  time  of  the  signal  out  of  the  antenna  increases  as  antenna  filtering  effects 

increase.  For  the  frozen-in  case,  C^t  +  Cyt  is  equal  to  unity  and  ta/To  =  'fQ-  This 
expression  gives  the  upper  limit  of  the  antenna-filtered  decorrelation  time. 
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UD 

Figure  11.  Filtered  decorrelation  time  for  a  uniformly-weighted  circular 
antenna  and  isotropic  scattering. 

The  corresponding  mean  Doppler  shift  is 

27caoCxt(Q-l)  ©o  4o  „om 

“A  =  ^  it'  (180) 

The  normalized  mean  Doppler  shift  due  to  antenna  pointing,  Xo°'A/Cxt .  is  plotted  in 
Figure  12  versus  the  ratio  Hq/D- 
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Figure  12.  Normalized  mean  Doppler  shift  due  to  antenna  pointing  for  a 
uniformly-weighted  circular  antenna  and  isotropic  scattering. 
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SECTION  4 

CHANNEL  SIMULATION 

A  statistical  channel  simulation  technique  is  described  in  this  section  that  allows 
realizations  of  the  impulse  response  functions  to  be  generated  at  the  outputs  of  multiple 
antennas  with  spatial  and  temporal  correlation  properties  given  by  the  GPSD  and  with 
Rayleigh  amplitude  statistics.  Realizations  generated  with  this  technique  represent  only 
the  diffractive  part  of  the  received  voltage,  and  they  are  valid  only  under  strong¬ 
scattering  conditions  where  the  GPSD  is  valid  and  where  Rayleigh  statistics  apply. 
Under  these  conditions,  however,  they  represent  a  solution  of  Maxwell's  equations  for 
propagation  of  RF  waves  through  randomly  structured  ionization. 

The  channel  simulation  technique  for  the  General  Model  that  is  described  here 
was  developed  by  Dr.  Leon  A.  Wittwer  of  the  Defense  Nuclear  Agency  (DNA),  and 
has  been  implemented  in  the  channel  model  Fortran  code  ACIRF  (Antenna/Channel 
/mpulse  Response  Function)  written  by  the  author  and  Dr.  Wittwer.  This  code  is 
available  to  qualified  users  through  DNA. 

The  basic  formalism  used  to  generate  statistical  realizations  of  the  channel 
impulse  response  function  without  antenna  effects  explicitly  included  was  first 
developed  by  Wittwer  [1980]  for  isotropic  irregularities.  Knepp  [1983]  extended  the 
technique  to  the  case  of  elongated  irregularities.  (The  elongated  case  corresponds  to  a 
90  degree  penetration  angle  and  to  an  infinite  axiai  ratio.)  The  channel  simulation 
technique  was  further  generalized  by  Dana  [1986]  to  include  the  effects  of  general 
anisotropic  scattering  and  multiple  antennas. 

4,1  IMPULSE  RESPONSE  FUNCTION. 

A  key  assumption  used  to  generate  realizations  of  the  channel  impulse  response 
function  is  that  the  channel  is  statistically  stationary  (in  the  wide  sense)  in  space, 
frequency,  and  time.  As  a  consequence,  the  impulse  response  function  is  delta 
correlated  in  angle,  delay,  and  Doppler  frequency.  This  allows  the  impulse  response 
function  to  be  generated  from  white  Gaussian  noise  in  the  angle,  delay,  and  Doppler 
frequency  domains,  and  then  Fourier  transformed  to  the  space  and  time  domains.  If 
necessary  for  a  particular  application,  the  Fourier  transform  from  delay  to  frequency 
may  be  performed  to  obtain  the  channel  transfer  function. 

The  impulse  response  function  at  the  output  of  an  aperture  antenna  located  at  p0 
and  pointed  in  direction  K0  is  given  by  taking  the  Fourier  transform  of  the  wave 
equation  solution  UA(po.w-t)  (Eqn.  121)  from  the  frequency  domain  to  the  delay 
domain  with  the  result: 

oo 

hA(p0A,t)  =  |  h(p,T,t)  A(p-p0)  exp  [iK0  (p-p0)j  d2p  .  (181) 

-OO 

This  equation  represents  the  spatial  convolution  of  the  aperture  weighting  function  and 
the  impulse  response  function  h(p,T,u  of  the  signal  incident  on  the  face  of  the  aperture. 
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The  channel  simulation  technique  depends  on  writing  h(p,T,t)  in  terms  of 
Fourier  transforms  from  angle  to  space  and  from  Doppler  frequency  to  time.  In  order 
to  show  that  Equation  181  can  be  written  in  this  form,  first  consider  writing  h(p,x,t)  in 
terms  of  its  Fourier  transforms: 

7d2K!  7dcoD  r  a 

h(p,x,t)  =  J  2  2  J  2k  exP  [>(Ki'P  -  coot)]  h(Ki,T,C0D)  •  (182) 

Equation  123  gives  the  aperture  weighting  function  A(p-p0)  in  terms  of  its  Fourier 
transform  and  the  antenna  voltage  pattern  g(Kj_).  Substituting  these  expressions  into 
Equation  181  and  doing  the  usual  change  in  the  order  of  integration  to  produce  a  delta 
function  results  in 


°°jl  oo 

r  U"Ki  r  dCOD  r  -»  a 

bA(po.Tt)  =  J  p  2  j  2k  exP  [‘(Ka  Po-wot)]  g(Ki-K0)  h(Ki,i,(OD)  . 

-OO  '  -oo 

(183) 

This  equation  is  used  to  generate  realizations  of  the  impulse  response  function  by 
first  generating  random  samples  of  ft(Ki,x,a)D)  and  then  performing  the  indicated 
Fourier  transforms.  Because  of  the  assumption  of  a  statistically  stationary  channel, 
(i(Ki,t,cod)  must  be  delta  correlated  in  angle,  delay,  and  Doppler  frequency: 

<fi(K^,T,ojD)ft*(K',T',mD))  =  S(K_l,t,g)D)  5(Ki-K’)  5(x-x')  5(cod-cod)  (184) 

The  first-order  amplitude  statistics  of  the  complex  quantity  h(p,x,t)  are  Rayleigh 
which  is  a  consequence  of  the  central  limit  theorem.  That  is,  h(p,x,t)  represents  the 
summation  of  many  scattered  waves  travelling  in  slightly  different  directions.  Thus  the 
two  orthogonal  components  of  h(p,x,t)  (either  the  in-phase  and  quadrature-phase 
components  or,  in  the  notation  used  in  this  report,  the  real  and  imaginary  parts)  are 
independent,  zero  mean,  normally-distributed  random  variables.  Consequently  the 
resulting  amplitude  is  Rayleigh  distributed,  and  the  resulting  phase  is  uniformly 
distributed.  Equation  181  indicates  that  hA(p,x,t)  is  the  summation  or  integration  of 
weighted  values  of  h(p,i,t)  and  is  therefore  also  Rayleigh  dis'ributed  (i.e.  the  sum  of 
normally-distributed  random  variables  is  itself  normally  distributed).  Strictly  speaking 
the  statistics  of  fi(K_,x,(i)D)  could  be  almost  anything  that  obeys  Equation  184,  and  the 
central  limit  theorem  could  be  invoked  to  argue  that  h(p,i,t)  and  hA(pTd)  are  zero- 
mean,  normally-distributed  complex  quantities.  Indeed,  multiple  phase  screen 
techniques  can  be  used  to  generate  Rayleigh-distributed  realizations  of  h(p,x,t)  starting 
with  just  random  phase  perturbations  of  an  electric  field.  However,  such  faith  in  the 
central  limit  theorem  is  not  necessary  if  fi(Ki,x,«D)  starts  out  as  a  zero-mean, 
normally-distributed  complex  quantity.  This  allows  many  fewer  points  to  be  used  in 
performing  discrete  Fourier  transforms  from  angle  to  space  than  would  otherwise  be 
required  to  guarantee  Rayleigh  statistics. 
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4.2  GENERATION  OF  REALIZATIONS. 


The  first  step  necessary  to  generate  realizations  of  the  impulse  response  function 
is  the  evaluation  of  the  GPSD  on  an  angular  Kx-Ky  grid.  Once  the  power  in  the 
angular  grid  cells  has  been  obtained,  these  quantities  are  used  to  construct  random 
samples  of  the  angular  spectrum  of  the  signal.  The  delta-function  relationship  between 
angle  and  delay  (Eqn.  105)  is  used  to  relate  an  annulus  in  the  Kx-Ky  grid  to  specific 
delay  bins,  and  the  translation  properties  of  the  GPSD  seen  in  Equation  105  are  used  to 
relate  specific  angles  to  Doppler  bins.  The  random  angular  spectrum  is  then  multiplied 
by  the  antenna  beam  pattern,  and  a  two-dimensional  discrete  Fourier  transform  (DFT) 
is  performed  to  the  antenna  phase  center  location.  A  final  Fourier  transform  from  the 
Doppler  frequency  domain  produces  the  impulse  response  function  as  a  function  of 
lime  and  delay.  Examples  of  such  realizations  are  presented  in  Section  5. 

Before  launching  into  a  detailed  discussion  of  the  channel  simulation  technique,  it 
is  useful  to  consider  the  consequences  of  the  shifting  property  of  the  angle-Doppler 
GPSD  (Eqns.  106  or  109). 

4.2.1  Computationally  Efficient  form  of  the  Impulse  Response 
Function. 


The  angle-Doppler  GPSD  in  Equation  106  contains  terms  of  the  form  Kx£x  - 
CxiT(jC0D  and  Ky  4y  -  CytT0C0D-  Thus,  except  for  the  leading  Gaussian  term  in  Equation 
106,  the  mean  signal  power  at  a  non-zero  Doppler  frequency  can  be  obtained  from  the 
GPSD  evaluated  at  zero  Doppler  frequency  and  shifted  in  angle  by  the  appropriate 
amount.  This  fact  is  used  to  reduce  the  computations  necessary  to  generate  realizations 
of  the  impulse  response  function. 

To  see  the  consequences  of  the  shifting  property,  consider  the  impulse  response 
function  given  by  Equation  183  which  is  the  basis  of  the  channel  simulation  technique. 
In  continuous  notation  and  using  the  delta-function  relationship  between  angle  and 
delay,  the  random  angle-Doppler-delay  spectrum  of  the  signal  may  be  written  as 


fi(K_L,T,C0D) 


i -  Ay(Kzx+K5)ix 

£n(Ki,C0d)  VSKD(K±,C0d)  if  T  =  -  — A ■-  . 


4(0coh 


(185) 


10 


otherwise 


where  sN(Kj_,0)d)  is  white  Gaussian  noise  with  unity  mean  power  and  with  the 
properties 

^N(Ki,o)D)qN(K',(ob)>  =  5(Ki-K')  5(coD  0)b)  (186a) 

<^N(Klto)D)^N(K',cob)>  =  0  •  (186b) 

Now  the  key  point  to  be  made  here  is  that  the  angle-Doppler  GPSD, 
SKD(Ex,Ky,0)D),  may  be  written  in  the  form 
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SKD(Kx,Ky,03D)  =  Sd(COd)  SKC(Kx-CxttoWD/ix,Ky-CytT0WD/iy) 


(187) 


where 


Sd(u>d)  =  V^T0exp 


_2 .  .2 

?pU>D 

4 


and  where 


SKC(Kx,Ky) 


7C&xfty 

Vi  -  Cxt  - 


(188) 


(189) 


x  exp 


KjAjO  -  Cv2,)  +  K2vt;(1  -  C>2,)4-2CxiCv,K,t>Kvllv 
4(1  -  Cxi  '  C2,) 


After  inserting  these  expressions  into  Equation  183,  the  impulse  response  function  is 

oo  oo  oo 

rdKx  r dKv  rdcon  r  -i 

hA(po*T,t)  =  J  ~2^  J  ~2^  \  ^-exp  [i(Kxx0+Kyyo-coDt)]  (190) 

-OO  -CO  -oo 

x  g(Kx-Kox,Ky-K0y)  ^N(KX,Ky,(0D) 

X  V  Sd(Wd)  SKC(Kx-CxtToWD/4x,Ky-CytT0COD/ly)" 

where  p()  =  (xo,y0).  Of  course  the  delta-function  relationship  between  angle  and  delay 
still  holds  although  it  is  not  shown  explicitly  in  this  equation.  The  problem  with  this 
expression  as  written  is  that  S«c  must  be  recomputed  for  each  new  Doppler  frequency, 
which  is  time  consuming. 

However  Skc  can  be  evaluated  once  at  zero  Doppler  frequency  and  shifted  for 
non-zero  frequencies.  In  a  digital  simulation  this  shift  is  most  efficiently  done  in 
discrete  steps.  Thus  let 


CxHoOd  -  mxAKx &x  +  £x &x 

CytTotoD  =:  myAKyHy  +  £yHy 


(191a) 

(191b) 


where  AKX  and  AKy  are  the  angular  grid  cell  sizes  that  will  be  used  to  numerically 
evaluate  Equation  190, 


mx  =  tnt 


f  CxiTpCOD 

L  IxAKx. 

I  CytT()tJ)D 
mv  =  int  i  „  rrr 

J  ^  ityAKy  J 


(192a) 

(192b) 
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and  ex  and  ey  are  the  residuals  after  the  discrete  shifts.  The  function  int  [■]  denotes  the 
integer  part  of  the  argument.  Now  define  a  shifted  angle-Doppler  GPSD,  denoted  Sks: 

SKS(Kx.Ky)  =  SKC(Kx"Cxt^o^D/^xiKy~CytTo£OD/£y)  .  (193) 

After  substituting  this  into  the  equation  for  the  impulse  response  function  and  changing 
angular  variables  to 

Kx  =  Kx-ex  (194a) 

Ky  =  Ky  -  £y  ,  (194b) 

Equation  190  becomes 

°°d  K '  °°d  K '  °°d 

hA(PoT4)  =  j  ~^r  J  J  ~|~exp  {i[(Kx+£x)x0+(I^+£y)yo-cODt]}  (195) 

-oo  -oo  -oo 

x  g(Kx+£x-K0x,Ky+£y-K0y)  ^N(Kx+£x,Ky+£y,G)D)  VSd(COd)  SKS(Kx>Ky)  • 

We  have  ignored  the  residual  shift  in  the  arguments  of  Sks  so  this  function  is  the 
result  of  a  discrete  shift  of  the  function  Skc-  Equation  195,  in  discrete  form,  is  used  to 
generate  the  impulse  response  functions  at  the  outputs  of  multiple  antennas. 

4.2.2  Discrete  Evaluation  of  the  GPSD. 

The  first  step  in  generating  the  impulse  response  function  at  the  output  of  an 
antenna  is  the  evaluation  of  the  GPSD  on  a  discrete  Kx-Ky-tOD  grid.  The  delta- 
function  relationship  between  angle  and  delay  is  used  to  relate  signal  components 
within  an  angular  annulus  to  a  particular  delay  bin.  Thus  at  this  point  it  is  not 
necessary  to  explicitly  include  delay,  and  the  GPSD  can  be  integrated  over  this 
variable.  In  order  to  assure  conservation  of  energy  (or,  more  correctly,  to  conserve 
signal  power),  the  GPSD  is  integrated  over  each  Kx-Ky-o)D  grid  cell,  and  that  power  is 
assigned  to  the  center  point  of  the  cell.  A  procedure  for  efficiently  performing  the 
three-dimensional  integral  is  described  in  this  subsection. 

The  two-dimensional  angular  and  DoppK°r  grids  are  defined  by  the  equations 

Kx  =  kxAKx  (-Nx/2<kx<  Nx/2-1)  (196a) 

Ky  =  kyAKy  (-Ny/2  <  ky  <  Ny/2-  1 )  (196b) 

cod  =  itidAcod  (-Nd/2  <  mo  <  Nd/2-1)  .  (196c) 

The  requirements  on  the  grid  cell  sizes  AKX,  AKy,  and  Acod  and  on  the  number  of  grid 
cells  Nx,  Ny,  and  Nd  will  be  discussed  later. 
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(kx+l/2)AKx  (ky+l/2)AKy  (mD+l/2)AcDD 
f  d K x  fdKv  f dcoo 

EA(kx,ky,mD)  =  I  2k  J  2  k  }  2k  sA(Kx,Ky,0)D) 

(kx-l/2)AKx  (ky-l/2)AKy  (mD-l/2)AcoD 

(197) 

Equation  197  is  completely  general,  but  it  implies  that  the  triple  integral  must  be 
computed  and  stored  separately  for  each  antenna  with  different  beamwidths  or  pointing 
angles.  If,  however,  it  is  assumed  that  the  antenna  beam  pattern  is  constant  across  a 
Kx-Ky-0L>D  grid  cell,  then  this  equation  can  be  approximated  by 

EA(kx,ky,mo)  =  G(kxAKx-K0x,kyAKy-K0y)  EKD(kx*ky,m[))  (198) 

where  EKD(kx.ky,mo)  is  the  incident  signal  power  in  the  Kx-Ky-C0D  grid  cell.  The 
accuracy  of  this  approximation  is  addressed  in  Appendix  C  where  it  is  shown  to 
conserve  energy  (or  power)  to  within  a  small  fraction  of  a  percent. 

In  order  to  produce  an  efficient  channel  model,  the  quantity  EKD(kx,ky,mD) 
must  be  readily  evaluated.  Therein  lies  a  major  problem  of  the  general  model. 
Straightforward  evaluation  of  Ekd  requires  a  closed-form  expression  for  the  triple 
integral: 

(kx+l/2)AKx  (ky+l/2)AKy  (mD+l/2)AtOD 
f  dKx  f  dKv  f  do)D  „ 

EKD(kx,ky,mD)  =  J  2K  J  2k  J  ~2k  sKD(Kx,Ky,tOD) 
(kx-l/2)AKx  (ky-l/2)AKy  (mD-l/2)AcoD 

(199) 

where  the  integrand  is  given  by  Equation  106.  It  is  clear  that  the  Kx-Ky  and  angle- 
Doppler  cross  terms  in  the  expression  for  SKD(Kx,Ky,0)D)  not  allow  a  simple  closed 
form  expression  for  Ekd  in  the  general  case,  although  such  expressions  can  be 
obtained  in  the  frozen-in  and  turbulent  limits. 

Two  “tricks”  are  used  to  efficiently  evaluate  Equation  199.  The  first  trick  is  to 
take  advantage  of  the  translational  properties  of  the  GPSD  described  in  the  previous 
subsection.  The  power  in  a  Kx-Ky-0)D  gnid  cell  is 

(mD+l/2)AcoD 

EKD(kx,ky,mD)  =  ^  2^  sd(o)d)  (200) 

(mD-l/2)Ao)D 

(kx+  1/2)AKX  (ky+l/2)AKy 
"  d  K  f  d  K 

X  ~  SKC(Kx'CxtT()(OD/(lx.Kv-Cvt'tou)D/ilv)  • 

I  J  A-iC  ' 

(kx-l/2)AKx  (ky-l/2)AKy 
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The  key  to  simplifying  this  expression  is  to  note  that  the  Doppler  grid  cell  size  is 
relatively  small  because  of  the  large  number  of  Doppler  samples  that  are  required  to 
produce  a  long  time  realization.  Thus  it  can  be  assumed  that  Skc  is  constant  over  a 
Doppler  cell,  and  Equation  200  reduces  to 

EKD(kx,ky,mD)  =  Eo(mD)  EKC(kx-mx,ky-my)  (201) 


where 


ED(mo)  =  2  ierfc 


(mD-l/2)ToAa)D~i 

2  J 


-  erfc 


(mp-i-l/2)T0A(0D 

2 


(202) 


The  quantity  EKc(kx-mx,ky-my)  is  the  power  in  a  shifted  Kx-Ky  grid  cell.  Because  of 
the  Kx-Ky  cross  terms  in  the  expression  for  Skc,  an  easily  evaluated  closed-form  result 
is  still  not  obtainable  for  Ekc. 


The  second  trick  used  in  the  channel  model  technique  is  to  note  that  a  rotation  by 
the  angle  0  (Eqn.  108)  in  the  Kx-Ky  plane  produces  an  orthogonal  form  of  the  GPSD 
which  does  not  contain  angular  cross  terms,  and  is  therefore  readily  integrated.  This 
orthogonalized  GPSD  is  given  by  Equation  109,  which  has  the  following  form  for  its 
angular  part: 


SKC(Kp,Kq)  = 


V(i  -  Cpt)( i  -  Cq2,)  expL  4<>  •  CP.)  4<>  -  CA)J 

The  signal  power  in  a  Kp-Kq  grid  cell,  with  indices  kp  and  kq  respectively,  is 
EKC(kP,kq)  =  Ep(kp)  Eq(kq) 

where 

(kp-f  1  /2)  AKp  S.  p 


KJtq 


(203) 


Ep(kp)  =  ^ \ erfc 


(kp-l/2)AKptp 

i.  2 VT~c|~. 


erfc 


2 V 1  -  Cp!, 


(204) 


A  similar  expression  holds  for  Eq(kq). 

Now,  EKC(kp,kq)  can  be  computed  on  a  fine  Kp-Kq  grid,  and  the  values  simply 
assigned  to  the  Kx-Ky  grid  cell  in  which  they  fall.  The  Kx-Ky  cell  indices  are 
computed  as  follows: 


kx 


int 


kpAKpCOSt)  -  kqAKqSint) 
AKx 


(205a) 


ky 


int 


kpAKpsinO  +  kqAKqCosO 
AKy 


(205b) 
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The  total  power  in  a  Kx-Ky  grid  cell  is  then  the  sum  of  all  EKC(kp-kq)  values  that  fall 
within  the  Kx-Ky  cell.  Roughly  ten  Kp-Kq  grid  cells  are  required  within  each  Kx-Ky 
cell  for  this  brute-force  procedure  to  work.  Thus  the  Kp-Kq  cell  sizes  are  determined 
by  the  expressions: 


AKp 


AKq 


0.1 


"  cos2t> 

sin20  ' 

L(AKx)2 

+  (AKy)2J 

0.1 _ 

sin20 

cos2d 

-(AKX)2 

+  (AKy)2J 

(206a) 


(206b) 


A  detailed  description  of  the  algorithms  used  to  compute  EKc(kx-irix.ky-my)  and 
to  shift  this  array  for  different  Doppler  frequencies  is  given  in  Appendix  D. 

4.2.3  Random  Realizations. 


The  next  step  in  the  channel  model  is  to  generate  a  random  realization  of  the 
angle-Doppler  spectrum  of  the  impulse  response  function  and  to  assign  the  spectral 
components  to  delay  bins  using  the  delta-function  relationship  between  angle  and  delay. 
Discrete  Fourier  transforms  are  then  performed  to  obtain  the  impulse  response 
function. 

4.2.3. 1  Assignment  of  Angular  Spectral  Components  to  Delay  Bins.  The  delta- 
function  relationship  between  angle  and  delay  (Eqn.  105)  in  the  diffraction-limited 
form  of  the  GPSD  is  used  to  assign  angular  spectral  components  to  discrete  delay  bins. 
This  form  of  the  GPSD  is  non-zero  only  when 


T 


Ay(KX  +  Ky)ll.X 


(207) 


A  straightforward  approach  to  assigning  angular  spectral  components  to  delay 
bins  is  to  compute  the  right-hand  side  of  Equation  207  at  the  center  of  each  Kx-Ky  grid 
cell  and  to  compute  the  index  of  the  delay  bin  using 


JL  =  Ay(K2x+K2v)l2x 
Ax  40)cohAx 


(208) 


where  At  is  the  sample  size  of  the  delay  bins.  The  problem  with  this  approach  is  that 
when  the  delay  sample  size  is  sufficiently  small,  the  number  of  angular  spectral 
components  that  fall  within  a  delay  bin  may  vary  substantially  from  one  delay  bin  to 
the  next  producing  ragged  statistics.  A  simple  solution  to  this  problem  is  to  randomly 
wiggle  the  angular  grid  cell  centers  before  applying  Equation  208.  This  spreads 
angular  spectral  components  more  or  less  uniformly  into  any  one  of  several  delay  bins, 
and  results  in  better  agreement  between  the  ensemble  signal  power  in  a  delay  bin  and 
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the  realization  power  in  that  bin.  The  randomly-wiggled  angular  grid  cell  centers  are 
computed  as 


Kx 


kx  +  ^u,x 


(209a) 


Ky 


ky  +  ^U.y 


(209b) 


where  ^u,x  and  ^u,y  are  independent,  uniformly  distributed  random  numbers  on  the 
interval  (0,1).  These  wiggled  cell  centers  are  then  used  to  compute  the  corresponding  j 
index  of  the  delay  grid. 

In  this  way  each  angular  spectral  component  is  assigned  to  a  unique  delay  bin. 
Because  the  angular  spectral  components  are  uncorrelated,  this  procedure  also 
guarantees  that  the  impulse  response  function  is  uncorrelated  from  one  delay  bin  to 
another. 


4. 2.3. 2  Discrete  Representation  of  Impulse  Response  Function.  The  discrete 
impulse  response  function  is  defined  on  a  discrete  time  and  delay  grid  defined  by  the 
equations: 

t  =  ntAt  (nt=l,  2,  -  ,  Nt)  (210a) 

t  =  JAt  .  0=0,  1,  -  ,Nt-1)  (210b) 

The  requirements  on  the  discrete  time  step  At,  the  number  to  time  steps  Nl(  the  delay 
sample  size  Ax,  and  the  number  of  delay  bins  NT  will  be  given  in  Section  4.3. 


Equation  195  gives  the  impulse  response  function  in  terms  of  Fourier  transforms 
from  the  random  angle-Doppler  spectral  components  to  the  time  domain  at  the  location 
of  the  phase  center  of  each  antenna.  In  discrete  form,  this  equation  is 


hA(jAx,ntAt) 


Nd/2-1  Nx/2-1  Ny/2-l 

8rc3 _ A  (Op  AKX  'V'1  AKy 

AKxAKyAcooAx  ZmU  2k  LmU  2k  2rt 

mD=-NTD/2  kx=-Nx/2  ky=-Ny/2 


x  exp  {i[(kxAKx+ex)x0+ (kyAKy+ey)y0  -  moAcoDntAt]}  (211) 


x  g(kxAKx+ex-K0X,kyAKy+ey-K0y)  ^N(kx,ky,mD) 
x  V Eo/mo)  EKS(kxAKx,kyAKy)  . 

The  normalization  factor  (87t3/AKxAKyAcODAx)  has  been  chosen  so  that  hA(jAx,ntAt)Ax 
represents  the  received  signal  during  the  delay  interval  jAx  to  (j+l)Ax.  As  the  Kx  and 
Ky  sums  are  performed,  Equation  208  is  used  to  assign  angular  spectral  components  to 
delay  bins.  Then  the  signal  components  in  each  delay  bin  are  Fourier  transformed 
from  the  Doppler  domain  to  the  time  domain. 
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The  quantity  ^\(kx,ky,mD)  is  a  complex,  zero-mean,  Gaussian  random  variable 
with  the  properties 

<'N(a,b,c)q\(a,p,Y)>  =  5a,a  5b.p  Sc,y  (212a) 

<'N(a,b,c))  =  0  (212b) 

(;N(a,b,c)CN(a,p,y))  =  0  (212c) 

where  8m,n  is  the  Kronecker  delta  symbol.  A  convenient  method  of  generating  the 
complex,  zero-mean,  Gaussian  random  numbers  is  to  use  the  following  equation 

=  y!-Po  In(qu.i)  exp  (2jci£u,2)  (213) 

where  2u,i  and  c.L’,2  are  independent  random  numbers  uniformly  distributed  on  the 
interval  [0,1 ),  and  P0  is  the  mean  power  of  the  random  samples  (P0  =  (c,_\cn)  =  1  for 
this  application). 

In  comparing  the  discrete  equation  for  the  impulse  response  function  with  its 
continuous-variable  analog  (Eqn.  195),  note  that  the  residual  shifts  ex  and  ey  may  be 
ignored  in  die  arguments  of  the  random  spectral  components,  £n.  because  shifted  white 
Gaussian  noise  is  still  white  Gaussian  noise. 

4. 2. 3. 3  Elimination  of  the  DC  Component.  As  written,  the  equation  for  the 
discrete  impulse  response  function  (Eqn.  211)  allows  a  Doppler  spectral  component 
with  zero-Doppler  frequency  (i.e.  the  mo  =  0  component).  This  component  will  in 
turn  result  in  a  DC  component  in  the  time  domain  realization,  which  is  undesirable, 
particularly  if  the  DC  component  is  large. 

A  simple  solution  to  this  problem  is  to  set  the  zero-Doppler  frequency 
component  of  c.\(kx,kx,mD)  to  zero.  Just  doing  this,  however,  results  in  reducing  the 
mean  power  in  the  impulse  response  function  by  the  zero-Doppler  frequency  power, 
ED(0). 

This  latter  problem  can  be  simply  solved  by  allowing  the  Doppler  frequency  bins 
adjacent  to  zero  Doppler  to  expand  in  size.  Thus  the  first  positive  Doppler  frequency 
bin  encompasses  frequencies  0  to  3Ag)d/2  and  has  power  EdO)  +  Ed(0)/2.  Similarly, 
the  first  negative  Doppler  frequency  bin  encompasses  frequencies  0  to  -3Ag)d/2  and 
has  power  EdCI  )  +  Ed(0)/2.  All  other  Doppler  frequency  bins  encompass  frequencies 
(mD- 1/2)Acod  to  (niD+l/2)Ao)D  and  have  power  Ed(itid). 

4.3  GRIDS. 

Angle  and  Doppler  frequency  grid  sizes  are  determined  by  requiring  that  the 
angular-Doppler  grid  encompass  a  large  fraction,  say  0.999,  of  the  power  in  the  GPSD 
of  the  signal.  The  0.001  error  must  then  be  allocated  between  the  angular  and  Doppler 
parts  of  the  three-dimensiona1  grid.  An  arbitrary,  but  intuitively  reasonable,  allocation 
is  to  divide  the  error  equally  between  the  angular  and  Doppler  frequency  parts  of  the 


6K 


grid  and  equally  between  the  two  angular  components.  Thus  the  Doppler  grid  limits 
are  determined  by  requiring  that  the  Doppler  grid  encompass  V0.999  of  the  Doppler 
frequency  power,  and  each  angular  grid  must  encompass  (0.999)'/4  of  the  angular 
power 


The  angular  and  Doppler  frequency  power  spectra  are  all  Gaussian.  Thus  each 
can  separately  be  written  in  the  form 


S(K) 


V  rr  exp 


(214) 


where  k  is  a  normalized  angle  or  Doppler  frequency  (i.e.  k  is  equal  to  KxUx  or  Kyfly 
or  TqCOd).  In  order  for  a  k  grid  to  encompass  a  fraction  £()  of  the  signal  power,  it  must 
extend  from  -Kmax  to  +Kmax  where 


Km  ax 

Co  =  |S(K)dK  .  (215) 

'Kmax 

This  equation  is  easily  solved  for  Kmax  in  terms  of  Co  with  the  result 

Kmax  =  2  erf-1  (Q  (216) 

where  erf-'  (•)  is  the  inverse  error  function.  If  t0  is  chosen  to  be  V 0.999  for  the 
Doppler  frequency  grid,  then 

KD.max  =  2x2.4612  =  4.9224  ,  (217a) 

and  if  Co  is  chosen  to  be  (0.999)'/4  for  either  angular  grid,  then 

KK.max  =  2x2.5895  =  5.1790  .  (217b) 

A  second  requirement  on  the  size  of  the  angular,  Doppler,  and  delay  grids  is  that 
they  be  defined  at  the  output  of  the  antennas,  thereby  eliminating  areas  of  tire  grids  that 
contribute  to  the  signal  power  incident  on  the  antennas  but  do  not  contribute  to  the 
power  of  the  output  signals  because  of  antenna  filtering.  Much  of  the  complexity  of 
the  algorithms  used  to  determine  grid  sizes  is  a  result  of  this  requirement,  but 
computing  grid  sizes  in  this  way  results  in  a  substantial  reduction  in  the  size  of  the 
grids,  and  therefore  in  computation  time,  when  antenna  filtering  effects  are  large. 

The  numbers  of  cells  in  the  angular  (Nx  and  Nv),  delay  (\x),  and  time  ( N , )  grids 
are  inputs  to  the  channel  simulation.  The  delay  sample  size  (At)  and  the  number  of 
samples  per  decorrelation  time  (N0)  arc  also  inputs.  From  these  quantities  and  the 
channel  and  antenna  parameters,  the  angular  (AKX  and  AKy),  Doppler  (Amo),  and  time 
(At)  grid  ceil  sizes  and  the  required  number  of  Doppler  samples  (No)  are  computed. 
Requirements  on  input  grid  parameters  and  consistency  checks  on  computed  grid 
parameters  are  described  in  the  next  subsections. 
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4.3.1  Angular  Grid. 


Angular  grid  sizes  are  determined  by  the  requirement  that  the  fraction  \  0.999  = 
0.9995  of  the  signal  power  after  antenna  filtering  be  contained  in  the  two-dimensional 
angular  grid.  First  consider  the  Kx  grid.  Because  of  the  symmetry  in  the  angular 
grid,  the  requirements  on  the  Ky  grid  can  then  be  obtained  by  analogy. 

The  Kx  power  spectrum  at  the  output  of  an  antenna  is 


dKv 


oo  oo 


d(0D 


-oo  -oo 


\  rcix 

VQy 


exp  ;  -  a uKp,u  -  orvK 


<)v 


whet . 


Ax  -  (Qx-DKox^x  +  OxyKoyly 
Ay  =  QxyKox  8-x  +  (Qy  1  )K0y  ®-y 
The  limits  of  the  Kx  grid  are  determined  by  requiring  that 


x-K 

XO.Ky-Kyo)  S(Kx,Ky,I.f0D) 

(218) 

Q5 

Kx*x  f  A  y  Q  x  y )  KXHX 

+  -Ai 

Qy 

4  "  +  ^Ax  Qy  j  ~i 

4Qy 

(219a) 

(219b) 

CoPa  = 


Kx.max 

dKx 

SA(KX)  ~2~ 


-K 


(220) 


x.max 


If  the  coefficient  of  the  linear  Kx  term  in  the  exponent  of  Equation  218  is  positive  then 
the  lower  limit  of  the  integral  in  Equation  220  can  be  replaced  by  -oo.  Conversely,  if 
the  coefficient  is  negative  then  the  upper  limit  can  be  replaced  by  Either  way,  the 
result  for  is 


r 

SO 


+  erf 


;  Q,>Kx  .max  H  x 

L  2Vq7 


lAxQy  -  AyQxyl 

2Q(,\rQy" 


(221) 


where  erf  (•)  is  the  error  function.  Setting  ^()  equal  to  (0.999) 1  and  solving  for 
Kx.max  gives  the  following  approximate  result: 


KK,max  lAxQy  -  AyQXyl 
+  4AxQ.,VQy 


(222a) 


Similarly,  the  limit  of  the  Kv  grid  is 
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(222b) 
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The  first  terms  in  the  expressions  for  KXimax  and  Ky.max  give  the  required 
extent  of  the  grid  when  the  antenna  is  pointed  along  the  line-of-sight.  The  second 
terms,  which  are  non-zero  only  when  the  antenna  is  pointed  away  from  the  line-of- 
sight,  give  the  amount  by  which  the  grid  must  be  extended  in  ordei  for  the  grid  to 
encompass  the  beam. 

Clearly,  Kx,max  and  Ky.max  depend  on  antenna  beamwidths  and  pointing  angles 
because  of  the  A  and  Q  factors.  If  there  are  multiple  antennas  then  Kx,max  and  Ky,max 
must  be  computed  for  each  antenna.  The  largest  values  are  then  used  to  determine  the 
boundaries  of  the  angular  grid. 


The  angular  grid  cell  sizes  can  now  be  computed  as 
2K, 


AK, 


AKV  = 


■  x.max 
NT 


2KV, 


max 


Nv 


(223a) 


(223b) 


where  the  number  of  grid  cells,  Nx  and  Nv,  are  inputs  to  the  channel  model. 

A  reasonable  minimum  value  for  the  number  of  Kx  or  Ky  grid  cells  is  32. 
However,  this  number  may  not  be  sufficient  if  there  are  multiple  antennas  with 
different  phase  center  locations.  Consider  the  two  antennas  with  the  largest 
separations,  dx  and  dy,  in  the  x-y  plane.  Because  the  impulse  response  functions  at  the 
outputs  of  the  two  antennas  are  generated  by  discrete  Fourier  transforms  from  the 
angular  domain  to  the  phase  center  locations  of  the  antenna,  the  unamb’r'  ious  distances, 
2n/AKx  and  27t/AKy,  of  the  DPI's  must  exceed  the  maximum  antenna  separations,  say 
h\  a  factor  of  2.  This  then  puts  upper  limits  on  AKX  and  AKy: 


AKX  < 


K 

dx 


(224a) 


AKV  <  f  .  (224b) 

Uy 

If  these  criteria  are  not  met,  then  \x  and/or  Ny  must  be  increased,  thereby  decreasing 
AKy  and/or  AKy,  until  they  are  met.  The  minimum  required  values  for  N-,  and  Ny  can 
then  he  written  as 


Nx 


max  32. 


s  x  ,m.i  x 
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( 225a) 
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(225b) 


4.3.2  Doppler  Frequency  and  Time  Grids. 

The  Doppler  frequency  grid  size  is  also  determined  by  the  requirement  that  99.9 
percent  of  the  signal  power  after  antenna  filtering  be  contained  in  the  Doppler 
frequency  grid.  Thus  the  antenna-filtered  Doppler  power  spectrum  is  required.  This 
function  is  most  easily  obtained  by  noting  that  the  temporal  coherence  function  has  the 
form 


T a ( t )  =  Pa  exp  - 


-i  +  i  m  a  t 

Ta 


(226) 


Recall  from  Section  3.2.3  that  coa  is  the  mean  Doppler  shift  due  to  antenna  pointing. 
The  Doppler  power  spectrum  is  then  given  by  the  equation 

7  -  r  ta  (wd+wa)2  " 

S.a(o)d)  =  I  F a ( t )  dt  =  v  k  ta  Pa  exp  -  — — q - -  :  .  (227) 


The  calculation  of  the  limits  on  the  Doppler  frequency  grid  is  exactly  analogous 
to  that  done  for  the  angular  grid.  The  limits  of  the  cod  grid  are  determined  by  setting 

equal  \  0.999  in  the  equation 


eOD.max 

'  do)D 

Ta  =  S  a  ( 0) d  )  ~2nT 

■tOD.max 


Solving  for 


1  „  “lA(O)D.max-ICi)Al) 

v-o  =  3  '  1  +  erf  - ^ 


(228) 


(229) 


gives  the  following  approximate  result  for  C0Q.max: 


tOD.max 


S.’  D.max 
TA 


+  I  CO  a! 


(230) 


The  first  term  in  this  expression  gives  the  required  maximum  Doppler  frequency 
when  the  antenna  is  pointed  along  the  line-of-sight,  and  the  second  term  is  the  result  ot 
die  mean  Doppler  shift  effect  of  antenna  pointing.  If  there  are  multiple  antennas,  then 
cod. max  must  be  calculated  for  each  antenna,  and  the  largest  value  used  to  determine  the 
Doppler  grid  size. 

At  this  point  the  required  number  of  Doppler  frequencies  Nd  is  still  unspecified. 
However,  because  a  fast  Fourier  transform  (FFT)  will  be  used  to  transform  from 
Doppler  frequency  to  time,  the  time  grid  requirements  may  be  used  to  derive  the 
Doppler  frequency  grid  cell  size  Acod  and  then  the  required  number  of  Doppler 
frequency  samples.  Nd. 
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Consider  the  requirements  on  the  time  samples.  Dana  (1982,  1 988 j  has  shown 
that  at  least  10  samples  per  decorrelation  time  are  required  to  accurately  reproduce  the 
temporal  statistics  of  Rayleigh  fading.  This  is  also  a  DNA  requirement  on  Rayleigh 
Fading  realizations  of  the  impulse  response  function  [Wittwcr  1980].  The  time  grid 
cell  size  is  then 


At  = 


"tA.min 


(231) 


where  N0  is  the  number  of  samples  per  decorrelation  time  (N0  must  be  greater  than  or 
equal  to  10),  and  TA.min  is  the  smallest  value  for  all  antennas  of  the  filtered 
decorrelation  time. 

In  addition,  DNA  requires  that  there  be  at  least  100  decorrelation  times  in  all 
realizations  of  the  antenna  output  impulse  response  function.  Thus, 


N.  > 


1 OOtA.max 

At 


(232) 


where  T..\,max  is  the  largest  value  for  all  antennas  of  the  Filtered  decorrelation  time.  If 
this  condition  is  not  met  then  the  number  of  time  samples  Nt  must  be  increased.  It  is 
also  necessary  that  Nt  be  equal  to  a  power  of  2  in  order  to  use  an  FFT.  The  minimum 
value  of  N’t  is  then  1024  in  order  to  meet  the  requirement  in  Equation  232  with  N0 
equal  to  10. 

Because  of  the  FFT  relationship  between  the  time  and  Doppler  frequency 
domains,  the  Doppler  frequency  grid  cell  size  is 


Amd  =  ■  (233) 

The  minimum  number  of  Doppler  frequency  samples  necessary  for  the  grid  to 
encompass  the  maximum  required  Doppler  frequency  is  then  given  by 


Nd  = 


2mp,max 

AMD 


‘^t^)D,max^A,min 

7TNo 


In  general,  N'd  will  be  smaller  than  Nt  implying  that  fewer  than  Nt  Doppler 
frequency  samples  are  required.  The  Doppler  frequency  arrays  may  then  be  zero- 
padded  to  N'i  samples  before  the  FFT  is  performed.  If,  however,  Np  is  greater  than 
N't,  the  implication  is  that  there  are  too  few  samples  per  decorrelation  time,  and  N() 
must  be  increased. 

Ihe  minimum  number  of  Doppler  frequency  samples  can  be  computed  from 
Equations  230  and  234.  If  all  antennas  are  pointed  along  the  iine-of-sight,  then  the 
minimum  value  of  Np  is 


’Sp.ruax-'si 


For  realizations  with  100  decorrelation  times  (Nt/N0  =  100),  the  minimum  number  of 
Doppler  samples  is  approximately  150.  If  antennas  are  pointed  away  from  the  line-of- 
sight,  the  required  number  of  Doppler  samples  will  increase  beyond  150. 

4.3.3  Delay  Grid. 

The  delay  sample  size  is  usually  chosen  on  the  basis  of  the  modulation  bandwidth 
of  the  transmitted  signal,  and  is  therefore  not  a  parameter  that  is  under  the  direct 
control  of  the  channel  simulation.  For  example,  in  a  phase-shift  keying  (PSK) 
application  there  must  be  at  least  two  delay  samples  per  channel  symbol  in  order  to 
accurately  simulate  the  transmitted  frequency  spectrum.  The  delay  sample  size  is  then 
generally  chosen  to  be  equal  to  one-half  the  channel  symbol  period.  In  a  frequency- 
shift  keying  (FSK)  application  with  frequency  hopping,  the  delay  sample  size  is  chosen 
so  that  the  unambiguous  frequency  bandwidth  of  the  impulse  response  function,  1/At, 
exceeds  the  hopping  bandwidth  by  some  comfortable  margin. 

The  number  of  delay  bins  is  an  input  to  the  channel  simulation.  The  requirement 
on  both  N't  and  At  is  that  the  realization  delay  grid  size  NxAt  encompass  at  least  97.5 
percent  of  the  delayed  signal  power  at  the  outputs  of  the  antennas.  A  somewhat  smaller 
percentage  is  used  to  define  the  limits  of  the  delay  grid  than  is  used  to  define  other  grid 
limits  because  of  the  slower  decay  of  signal  power  with  delay  than  with  angle  and 
Doppler  frequency  (exponential  decay  with  delay  versus  Gaussian  decay  with  angle  and 
Doppler  frequency). 

'Ihe  ensemble  signal  power  in  the  delay  bins  is  given  by  the  integral 
(j+ 1 )  At 

Pj  =  J  SA(t)  dT  (236) 

jAt 

where  Sa(t)  is  the  delay  power  spectral  density  at  the  output  of  an  antenna.  The 
general  expression  for  Sa(G  is  quite  complicated  when  the  antenna  is  pointed  away 
from  the  line-of-sight,  and  will  not  be  given  here.  The  reader  is  referred  to  Frasier 
1 1988]  for  details  on  Sa(t). 

'I’he  total  signal  power  in  the  delay  grid, 

Nt-1 

PGt  =  I  Pj  ,  (237) 

j=0 


must  be  greater  than  or  equal  to  0.975Pa-  If  not,  then  either  Nx  or  At  or  both  must  be 
increased. 
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SECTION  5 

MATCHED  FILTER  EXAMPLES 

This  section  presents  examples  of  the  received  voltage  out  of  a  filter  matched  to 
a  transmitted  square  pulse.  These  examples  are  intended  to  illustrate  the  effects  of 
frequency  selectivity  and  antenna  filtering  on  a  transionospheric  communications  link 
and  to  illustrate  the  differences  in  the  structure  of  the  received  signal  depending  on 
whether  the  frozen-in,  turbulent,  or  general  models  are  used  to  generate  the  impulse 
response  function  realizations.  The  following  calculation  also  illustrates  how  the 
received  voltage  can  be  constructed  from  the  impulse  response  function  realizations  in 
a  digital  link  simulation.  Additional  examples  for  specific  system  applications  may  be 
found  in  Bogusch,  et.  al.  [1981  ]  and  in  Bogusch,  Guigliano,  and  Knepp  [19831. 

5.1  MATCHED  FILTER  OUTPUT  SIGNAL. 

The  output  of  a  matched  filter  can  be  constructed  by  convolving  the  impulse 
response  function  of  the  channel  and  antenna  with  the  combined  impulse  response 
function  of  the  transmitter  and  receiver.  A  second  approach  is  to  construct  the 
combined  frequency  response  of  the  transmitter,  channel,  antenna,  and  receiver  and 
then  to  Fourier  transform  that  result  to  obtain  the  matched-filter  output.  This  latter 
approach  is  used  for  the  examples  presented  in  this  section. 

The  starting  point  is  to  calculate  the  channel/antenna  transfer  function  which  is 
the  Fourier  transform  of  the  impulse  response  function: 


H(o),t)  =  j  h(x,t)  exp  (-ion)  di  .  (238) 

0 

This  function  represents  the  response  of  the  channel  and  antenna  at  time  t  to  a 
transmitted  sinew'ave  with  radian  frequency  0). 

For  a  transmitted  square  pulse  with  a  chip  duration  Tc,  the  voltage  out  of  the 
matched  filter  at  time  t  can  be  written  as 


i 


U(T,t) 


J  M(co)  H(o),t)  exp  (io)X) 

-oo 


do) 

2k 


(239) 


where  i  is  the  time  delay  of  the  matched  filter  relative  to  the  nominal  time-of-arrival 
(i.e.  the  time-of-arrival  under  benign  propagation  conditions).  The  combined 
spectrum  of  the  transmitted  square  pulse  and  the  receiver  matched  filter  is 


M(0)i 


sin:(o)Tc/2 ) 
(wTt/2): 


(240) 


The  impulse  response  function  is  generated  with  delay  samples  ot  size  At.  so 
the  discrete  channel/antenna  transfer  function  has  an  unambiguous  frequency  response 
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of  2ti/At  radians.  If  this  bandwidth  is  divided  into  Nf  frequency  samples,  then  the 
discrete  channel/antenna  transfer  function  at  time  ntAt  is 


N'x-1 

H(kFAco,ntAt)  =  £  h(jAx,ntAt)Ax  exp  (-ijAxkFAo)]  (241) 

]=0 

where  Am  =  2rc/(N'FAx).  Recall  that  the  normalization  of  the  impulse  response 
function  is  such  that  the  factor  At  following  h(jAx,ntAt)  must  be  included.  The  range 
of  the  index  kp  in  this  equation  is  from  -Nf/2  to  Nf/2-1  representing  a  range  of 
frequencies  from  -NpAm/2  to  (Ntf-1)Aco/2.  Of  course  the  number  of  frequency 
samples  must  be  at  least  as  large  as  the  number  of  delay  samples  in  order  for  the 
transfer  function  to  preserve  the  information  contained  in  the  impulse  response 
function.  However,  if  the  number  of  delay  samples  is  chosen  to  be  the  minimum 
number  required,  then  it  may  be  necessary  to  select  the  number  of  frequency  samples 
to  be  larger  than  the  number  of  delay  samples  in  order  to  minimize  aliasing  (in  delay) 
of  the  matched-filter  output. 

The  output  voltage  of  the  matched  filter,  as  a  function  of  time  and  relative  delay, 
is  then  given  by 


Ao)Tc 

u(x,ntAt)  =  ~2n~ 


N'f/2-1 

1 A ! 

kp-Np/2 


sin2(kFAcoTc/2)  j 


(kFAo)Tc/2)2 


2  H(kFAo),ntAt)  exp  (ikpAtox)  . 


(242) 


If  the  delay  samples  size  Ax  of  the  realization  of  the  impulse  response  function  >s 
chosen  to  be  Tc/2,  then  AorTc/2  =  2tx/Nf  and  u(x,n[At)  represents  a  signal  that  is  band- 
limited  to  the  frequency  range  -1/TC  to  +I/TC.  Note  that  the  matched-filter  output 
uix,t)  is  unambiguous  in  delay  over  the  interval  from  0  to  (Nf-1  )At  compared  to  the 
delay  interval  of  0  to  ( Nx- 1  )Ax  for  the  original  realization. 


In  the  examples  that  follow,  the  chip  rate  is  set  at  1  MHz,  and  the  random 
realizations  of  the  impulse  response  function  are  generated  with  a  delay  sample  size  of 
TJ2  (Tc  =  1/RC).  However,  the  frequency  selective  effects  depend  only  on  the  ratio  of 
the  frequency  selective  bandwidth  to  the  chip  rate  f,/Rc.  For  the  antenna  examples,  a 
uniformly-weighted  circular  antenna  and  isotropic  scattering  are  assumed.  Antenna 
effects  then  depend  only  on  the  ratio  of  the  antenna  diameter  D  to  the  decorrelation 
distance  and  the  antenna  pointing  angle. 

5.2  FRKQl  KM  V  SKLKCTIVK  EFFECTS. 

In  a  high  data  rate  communications  link,  the  major  effect  of  frequency  selective 
fading  is  intersymbol  interference.  Kven  relatively  small  amounts  of  delay  spread  can 
catastrophically  degrade  demodulation  performance  in  such  a  link  using  conventional 
matched-filter  detection  techniques. 
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Figure  13  shows  examples  of  the  matched-filter  output  amplitude  for  three  levels 
of  frequency  selective  propagation  disturbances,  characterized  by  the  ratio  of  the 
frequency  selective  bandwidth  f0  to  the  chip  rate  Rc.  The  impulse  response  functions 
were  generated  using  the  frozen-in  model  (Cxt  =  1  -0  and  Cyt  =  0)  and  a  small  antenna 
(i.e.  D  «  £0).  Each  frame  in  the  figure  provides  a  three-dimensional  picture  of  the 
matched-filter  output  amplitude  for  a  single  transmitted  pulse  as  a  function  of  time 
delay  (abscissa)  and  time  (scale  directed  into  the  figure).  The  total  duration  of  each  of 
the  frames  is  10  decorrelation  times. 

In  the  top  frame  the  frequency  selective  bandwidth  is  equal  to  the  chip  rate  and 
only  a  small  amount  of  distortion  is  evident  in  the  waveform  (which  is  slightly  rounded 
due  to  band  limiting  at  the  first  nulls  of  the  signal  spectrum).  The  effect  of  fading  can 
be  seen  in  this  frame  as  the  peak  amplitude  rises  and  falls  with  time.  Some  minor 
distortion  of  the  output  amplitude  is  seen  but  for  the  most  part  the  signal  is  contained 
within  the  period  of  one  chip.  This  channel  is  nearly-flat  fading  which  means  that  all 
frequency  components  within  the  signal  bandwidth  propagate  essentially  the  same  way 
through  the  disturbed  ionosphere.  There  is  very  little  time  delay  spread  beyond  one 
chip  in  the  matched-filter  output. 

The  middle  frame  in  Figure  13  shows  the  matched-filter  output  amplitude  for 
the  case  were  fo/Rc  is  equal  to  0.2.  For  this  smaller  value  of  the  frequency  selective 
bandwidth,  more  of  the  signal  energy  is  arriving  with  delays  of  more  than  a  chip,  and 
there  are  multiple  distinct  peaks  in  the  matched-filter  output  amplitude.  It  is  these 
structures  that  can  cause  delay  tracking  algorithms  to  lose  lock  and  that  cause 
intersymbol  interference  which  can  degrade  demodulation  performance. 

The  bottom  frame  shows  a  highly  disturbed  case  where  f0  is  a  tenth  of  the  chip 
rate.  This  causes  signal  energy  to  be  spread  over  approximately  eight  chip  periods. 
When  a  contiguous  set  of  pulses  is  transmitted,  the  delay  spread  of  the  received  signal 
results  in  the  simultaneous  reception  of  information  from  about  eight  previous  chips 
w'hich  can  produce  severe  intersymbol  interference. 

An  effect  due  to  the  frozen-in  model  that  is  evident  in  Figure  13  is  that  the  signal 
arriving  at  long  delays  varies  more  rapidly  in  time  than  the  signal  arriving  at  shorter 
delays.  A  comparison  of  the  matched-filter  output  amplitude  generated  with  the 
frozen-in,  general,  and  turbulent  models  is  shown  in  Figure  14  for  the  case  where  f(/Rc 
is  equal  to  0.1.  The  top  frame  in  this  figure  for  the  frozen-in  model  is  just  a 
reproduction  of  the  bottom  frame  in  Figure  13.  Again  10  decorrelation  times  of  the 
signal  are  plotted.  The  middle  frame  is  a  general  model  realization  (Cxt  =  0.9  and  Cy( 
=  0),  and  the  bottom  frame  is  for  the  turbulent  model  (Cxt  =  Cyt  =  0).  The  difference 
between  the  top  and  bottom  frames  is  that  the  turbulent  model  amplitude  has  the  same 
fading  rate  at  all  delays.  It  can  be  seen  that  the  general  model  realization  falls 
somewhere  between  these  two  limiting  cases. 
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Figure  13.  Effects  of  frequency  selective  fading. 
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Figure  14.  Comparison  of  matched  filter  output  amplitude  for  the  frozen- 
in,  general,  and  turbulent  models. 
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5.3  SPATIALLY  SELECTIVE  EFFECTS. 

Spatially  selective  effects  are  important  for  high  data  rate  communications  links 
that  rely  on  large  antennas  to  achieve  sufficient  signal-to-noise  ratios  for  low  error  rate 
data  demodulation.  Scattering  loss  of  an  antenna  is  a  function  of  the  size  of  the  antenna 
D  relative  to  the  decorrelation  distance. 

When  40  is  greater  than  D,  the  electric  field  is  highly  correlated  across  the  face 
of  the  antenna  and  the  full  gain  of  the  antenna  is  realized.  However,  the  antenna  may 
be  located  at  a  position  where  the  incident  power  is  in  a  deep  fade.  The  solution  to  this 
problem  is  to  have  multiple  antennas  physically  separated  by  a  distance  larger  than  the 
maximum  decorrelation  distance.  The  probability  of  having  all  antennas 
simultaneously  experience  deep  fades  in  the  received  power  is  then  substantially 
reduced. 

The  problem  of  spatial  selectivity  occurs  when  &0  is  less  than  D  and  the  electric 
field  is  decorrelated  across  the  face  of  the  antenna.  In  this  case,  the  induced  voltages  in 
the  antenna  add  noncoherently  due  to  the  random  phase  variations  in  the  electric  field, 
and  a  loss  of  signal  power,  or  equivalently  of  antenna  gain,  is  the  result.  From  another 
perspective,  this  loss  occurs  when  the  angular  scattering  process  responsible  for 
amplitude  and  phase  scintillation  and  frequency  selective  effects  also  causes  some  of  the 
transmitted  signal  energy  to  be  scattered  out  of  the  antenna  beam. 

Figure  15  shows  examples  of  the  matched-filter  output  amplitude  for  three  levels 
of  spatially  selective  propagation  disturbance,  characterized  by  the  ratio  of  the  antenna 
size  to  the  decorrelation  distance.  The  ratio  of  the  frequency  selective  bandwidth  to  the 
chip  rate  rate  is  0.1,  and  the  antenna  is  pointing  along  the  line-of-sight.  The  top  frame 
is  for  the  case  where  !0  is  much  greater  than  D,  and  is  just  a  reproduction  of  the 
bottom  frame  of  Figure  13.  The  middle  frame  is  for  a  i0/D  ratio  of  0.5  where  the 
scattering  loss  is  3.1  dB.  The  effect  of  the  antenna  is  to  preferentially  attenuate  the 
signal  energy  arriving  at  large  angles  and  also  at  large  delays  and  thereby  to  reduce  the 
delay  spread  of  the  output  signal.  In  the  bottom  frame  where  St(/D  is  equal  to  0.2,  the 
output  signal  is  almost  flat  with  very  little  delay  spread  distortion  of  the  matched-filter 
output.  Although  this  substantially  reduces  the  effects  of  frequency  selective  fading, 
the  cost  is  an  8.8-dB  reduction  in  the  average  signal  power. 

Finally,  Figure  16  shows  examples  of  the  matched-filter  output  amplitude  for 
three  values  of  the  pointing  angle  ©0.  The  ratio  of  the  frequency  selective  bandwidth 
to  the  chip  rate  rate  is  0.1,  and  the  ratio  of  decorrelation  distance  to  antenna  diameter 
is  0.5.  The  top  frame  for  a  pointing  angle  of  zero  is  just  a  reproduction  of  the  middle 
frame  of  Figure  15.  The  average  scattering  loss  for  this  case  is  3.1  dB.  The  bottom 
two  frames  show  the  matched-filter  output  amplitude  for  pointing  angles  of  one-half 
beamwidth  (0o=9o/2)  with  a  scattering  loss  of  4.6  dB  and  one  beamwidth  (0O=0O)  with 
a  scattering  loss  of  9.2  dB. 
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Figure  15 
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Effects  of  spatially  selective  fading. 
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Although  the  average  scattering  losses  are  about  the  same,  the  bottom  frame  of 
Figure  16  (!</D  =  0.5  and  ©O=0O)  and  the  bottom  frame  of  Figure  15  ( i()/D  =  0.2  and 
@o=0)  are  qualitatively  quite  different.  For  the  case  with  the  pointing  angle  equal  to  a 
beamwidth,  the  received  power  is  much  more  spread  out  in  delay  compared  to  the  case 
with  zero  pointing  angle  where  the  signal  energy  is  concentrated  near  zero  delay.  This 
is  due  to  the  fact  that  the  antenna  pointed  away  from  the  line-of-sight  has  relatively 
higher  gain  at  larg .  angles  and  long  delays  and  relatively  lower  gain  at  small  angles 
and  short  delays  than  does  an  antenna  pointed  along  the  line-of-sight.  Thus  for  an 
antenna  pointed  away  from  the  line-of-sight,  increased  scattering  loss  does  not 
necessarily  result  in  reduced  frequency  selective  effects. 
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Figure  16.  Effects  of  beam  pointing. 
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APPENDIX  A 

PHASE  VARIANCE  DUE  TO  ELECTRON  DENSITY  FLUCTUATIONS 


A  relationship  between  the  phase  variance  imparted  on  the  wave  as  it  propagates 
through  the  ionization  layer  and  the  electron  density  fluctuations  will  be  derived  in  this 
appendix.  This  relationship  is  given  by  Equation  51  which  was  derived  using  the 
Markov  approximation.  However,  it  will  be  shown  here  that  the  relationship  requires 
only  that  the  layer  thickness  be  large  compared  to  the  decorrelation  distance  of  the 
electron  density  fluctuations  along  the  line-of-sight. 

It  was  showm  in  Section  2  that  the  total  phase  change  of  the  wave  as  it  propagates 
through  the  ionization  layer  is 


<t>  =  rcX(ne>  Js(p,z,t)dz 


(243) 


The  autocorrelation  function  of  the  phase  fluctuations  is  then 

L  L 

<<J)(p.t)<j>(p',t’))  =  (reX(nc))2  J  dz  J  dz’  (£(p,z,t)^(p',z',t')>  (244) 


where  L  is  the  thickness  of  the  scattering  layer.  For  spatially  and  temporally  stationary 
random  electron  density  fluctuations,  the  expectation  must  be  a  function  of  the 
differences  8p=p-p\  8z=z-z’  and  8t=t-t'  only.  Denoting  the  autocorrelation  of  £(p,z,t) 
by  B^(8p,8z,St)  and  the  autocorrelation  of  <J>(p,t)  by  B<j,(8p,8t),  Equation  244  becomes 

L  L 

B<j)(8p,8t)  =  (rcX(nc))2  j  dz  j  dz'  B^(5p,8z,8t)  .  (245) 

0  0 

This  double  integral  may  be  reduced  to  a  single  integral  by  changing  the  order  of 
integration  with  the  result 


L 


B$(8p,8t)  =  (reX{ne»2  L  j  dz 

-L 


( 

1 


B^(Sp,Sz,8t)  . 


(246) 


If  the  correlation  distance  of  B^(8p,8z,8t)  along  the  z  direction  is  small 
compared  to  L,  then  B?(8p,8z,8t)  will  become  small  before  Izl/L  approaches  unity  in 
the  integral,  and  the  Izl/L  term  may  be  ignored.  The  limits  of  the  integral  may  then  be 
set  to  +  °°  and  the  integral  reduces  to 


00 

Bo(6p,8t)  =  (rcX(nc))2  L  J  dz  Bc,(8p,8z,8t) 

-00 


(247) 
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The  remaining  integral  is  denoted  by  A^(5p,5t)  so  the  autocorrelation  of  the  phase 
fluctuations  is 

B0(5p,5t)  =  (rcX(nc))2  L  A^(Sp,5t)  ,  (248) 

which  is  the  same  as  the  final  expression  in  Equation  51. 
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APPENDIX  B 

CHANNEL  PARAMETERS  FOR  K  4  ELECTRON  DENSITY 

FLUCTUATIONS 


The  expansion  coefficients  A0  and  A2  are  calculated  in  this  appendix  using  the 
quadratic  approximation  of  the  structure  function  A^(p)  of  the  electron  density 
fluctuations  and  using  the  delta-layer  approximation.  From  these  coefficients,  the 
phase  variance,  decorrelation  distance,  and  the  coherence  bandwidth  of  the  signal 
incident  on  the  plane  of  the  receiver  are  written  in  terms  of  physical  parameters. 
However,  these  channel  per  imeters  are  computed  from  a  disturbed  ionosphere  model 
using  the  more  general  formalism  of  Wittwer  [1979,1980]  which  accounts  for  the  finite 
thickness  of  the  scattering  region  and  other  complicating  effects.  The  purpose  of  this 
appendix  is  only  to  illustrate  the  dependence  of  the  channel  parameters  on  geometrical 
and  electron  densitv  fluctuation  parameters. 

A  power-law'  form  of  the  PSD  for  the  three-dimensional  electron  density 
fluctuations  is  assumed: 


8tt3/2r(n)LrLsLt  (An2) 
T(n-3/2)  <ne>2  (1  +  K  L  K)n 


(249) 


where 


L 


"L2X  0  0  " 

0  Ly  LyZ 

_  0  Lyz  L2  _ 


(250) 


The  scales  Lr ,  Ls ,  Lt  and  Lx  ,  Ly  ,  Lz  ,  Lyz  are  defined  in  Section  2.6  in  terms  of  the 
penetration  angle  O  and  the  axial  ratio  q.  For  K  L  K  »  1,  S£  is  proportional  to  K'2n. 

Thus  a  K  4  PSD  for  the  three-dimensional  electron  density  fluctuations  corresponds  to 
the  n  =  2  case. 


The  structure  function  A^(p)  under  the  delta-layer  approximation  is  given  by 
Equation  71  which  is  reproduced  here: 


7  d2Ki 

A^(p)  =  J  exp  (iKi-p)  S^(Ki,Kz=0)  • 


(71) 


Using  Equations  249  and  250  and  performing  the  angular  integral.  Equation  71 
reduces  to 
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CO 


4VrtT(n)  qL~x(An 2)  Jo(K  V  P'Li'P 
A^(p)  =  F(n-3/2)  Ly  (nc>2  J  (1  +  K2)n 


K  dK 


(251) 


Vrc(n-l)  qLx  (An;;) 
2n-3F(n-3/2)  Ly  <ne) 


5(p  L-|,  p,«n'l>'2KI1M(Vp  L)  P  ) 


where  J0  is  the  Bessel  function  of  order  0,  Kn-!  is  the  modified  Bessel  function  of  order 
n-1,  Lj_  is  a  2-2  matrix  containing  the  x-y  components  of  L,  and  L'|  is  the  inverse 

matrix  of  Li. 

For  all  values  of  n  except  n  equal  to  2,  A^(p)  can  be  expanded  in  a  power  series 
of  the  form 


T  x2  y2  1  2 

A^(p)  -  A()  i  1  -  A  2  I  ,2  +  ,2 

[  L  Lx  Ly 


(252) 


where  m  =  min  (2,2n-2).  For  the  n  equal  2  case,  A2  does  not  exist  unless  an  inner  scale 
5.,  is  imposed.  This  is  accomplished  by  truncating  the  integral  over  K  in  Equation  251 
at  a  cutoff  K  -  Lx/!j.  For  values  of  n  greater  than  2,  the  J0  Bessel  function  in  the 
integrand  of  Equation  252  can  be  expanded,  and  the  resulting  series  can  be  integrated 
term-by-term.  The  first  twu  terms  of  the  expansion  give 


An  = 


2\^Ttr(n)  qL2  (An2) 
F(n-3/2)  Ly  <ne>2 


2  ll  -n 


1  A 
+  If 


(253) 


A2  = 


2  in- 1 


1  +  (n-1) 


2  in- 1 


4(n-2)|l  +  ^fJ  -1 


(254) 


In  the  limit  that  n  equals  2,  the  A2  coefficient  becomes 


A->  = 


In  ;  1  + 


(255) 


This  expression  can  be  further  reduced  in  the  limit  that  the  inner  scale  is  much  smaller 
than  the  outer  scale  (i.e.  ij  «  Lx)  to  give 
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(256) 


In  (Lx/4,) 

A2  -  j - 

Now  the  phase  variance,  decorrelation  distances,  and  coherence  bandwidth  can 
be  written  in  terms  of  geometrical  parameters  and  electron  density  fluctuation 
parameters.  Using  Equation  52,  the  phase  variance  due  to  structured  ionization  is 


2qLtL 


Lx 


rc^2  (An;:) 


(257) 


where  L$  is  the  thickness  of  the  delta  layer.  The  decorrelation  distances,  decorrelation 
time,  and  coherence  bandwidth  are  defined  in  Equations  79,  80  and  81  respectively. 
Using  these  definitions, 


and 


4*  = 


iy  = 


*o  = 


-i  i 


ln(Lx/4i) 

2 

ln(Lx/i  i) 


2  Lx  Zt  +  zr 
CJ<()  zt 

2  Ly  ?t±_zjr 
Cty  zt 


i  i 


ln(Lx/4i) 


2  lo 

0(j> 


wcoh  -  -> 


AyC0g 


Lx  zt  +  zr 


(258a) 

(258b) 

(259) 


c  ln(Lx/4i)  oj  ztzr 


(260) 


These  equations  are  only  valid  for  the  delta-layer  approximation,  for  the 
quadratic  phase  structure  approximation,  and  for  a  K'4  three-dimensional  electron 
density  fluctuation  PSD.  It  can  be  seen  from  the  equations  that  the  decorrelation 
distances  and  coherence  bandwidth  are  only  weakly  dependent  on  the  inner  scale. 
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APPENDIX  C 

ACCURACY  OF  ANGULAR  INTEGRATION  TECHNIQUES 


C.l  INTRODUCTION. 

The  mean  signal  power  in  a  Kx-Kv  grid  cell  at  the  output  of  an  antenna  is 
<kx+l/2)AKx  (ky-l/2)AKy 

H.\(kx,kx)  =  q_/  q:T  G( ,KX.KV)  Sk(Kx.Kv )  •  (261) 

ZJZ 

(kx-l/2)AKx  (ky-l/2)AKy 


Equation  261  is  quite  general,  but  it  requires  that  the  power  in  each  angular  bin  be 
calculated  and  stored  for  each  antenna.  This  latter  requirement  results  in 
unacceptablely  large  arrays.  However,  Equation  261  can  be  approximated  by  assuming 
that  the  antenna  beam  pattern  varies  slowly  over  the  grid  cells  so  G(KX.KV)  may  be 
pulled  out  of  the  integral.  ITiere  are  several  ways  that  this  can  be  done.  The  purpose 
of  this  appendix  is  to  calculate  the  accuracy  of  a  few  approximations  to  Equation  261. 

In  order  to  limit  the  scope  of  this  calculation,  isotropic  scattering  and  a 
uniformly-weighted  circular  antenna  will  be  assumed.  Without  further  loss  of 
generality,  it  can  then  be  assumed  that  the  antenna  is  pointed  away  from  the  line-of- 
sight  in  the  x-direction,  or  equivalently,  that  the  pointing  azimuth  is  zero.  For  this 
case,  the  angular  part  of  the  GPSD  is 


Sk(Kx,Kv)  =  exp 


K: )•?  1 


(262) 


and  the  antenna  beam  pattern  is 


G(KX,KV)  =  exp 


(Q-l  UK.-Kori,-,  (Q-l  )KUr; 
4  '  4 


where  Q-l  is  proportional  to  the  square  of  the  ratio  of  the  antenna  diameter  D  to  the 
decorrelation  distance: 


0 


4in2  D2 
(1.028997t):  if, 


(264) 


For  this  isotropic  scattering  and  antenna  case,  the  power  at  the  output  of  the 
antenna  reduces  to 


where  P>0  is  the  pointing  direction  elevation  angle  and  0()  is  the  Mam  width  of  the 
antenna. 
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The  antenna-filtered  decorrelation  distance,  which  is  the  same  in  both  the  and 
y-directions.  is 


l.\  =  .  (266) 

When  the  antenna  is  pointed  in  the  Kx  direction,  the  Kx  angular  grid  size  is  then  given 
by  die  expression 


2  (Q'ljKjj  ^  K  .max 

N'x  Q  +  40\Q 


( 267a) 


and  die  Ky  angular  grid  size  is 


AKy 


2  '  JyKjiiax  ; 

N>  :  a0\ q  j 


(267b) 


where  KK.max  (Eqn.  217b)  is  determined  by  the  condition  that  99.9  percent  of  the 
signal  energy  be  in  the  Kx-Ky  grid. 

Ilie  exact  expression  for  the  received  power  will  be  compared  to  the  total  power 
in  the  grid, 


N'x/2-1  Ny/2-1 

PGK  =  X  X  EA(kx,ky)  , 

kx=-N*/2  ky=-Ny/2 

to  conipuie  an  error  in  the  total  power 


PE 


p,\  -  Pgk 

!  Pa 


(268) 


(269) 


for  each  of  the  algorithms  used  to  evaluate  Equation  261 . 

C.2  ALGORITHMS. 

The  first  approximation  to  the  exact  result  is  just  Equation  261.  This  expression 
will  result  in  some  error  in  the  total  power  because  of  round-off  errors  in  the 
summation  of  the  contributions  from  each  grid  cell  and  because  of  the  finite  size  of  the 
angular  grid.  Additional  error  results  when  the  beam  is  pointed  near  the  edge  of  the 
grid  so  part  of  the  beam  is  pointed  out  of  the  grid.  The  magnitude  of  this  error  will  be 
apparent  in  the  results  of  this  analysis.  With  the  assumptions  of  isotropic  scattering  and 
an  isotropic  Gaussian  beam  pattern,  the  integrals  indicated  in  Equation  261  can  be 
obtained  in  closed  form  with  the  result 
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A 


The  advantage  of  this  approximation  is  that  ihe  error  function  terms  depend  oniy  on 
the  environment  so  they  can  be  done  once  and  used  for  all  antennas  thereby  reducing 
the  required  processing  time  and  array  sizes. 

A  third  approximation  to  Equation  2b  1  is  similar  to  the  previous  approximation. 
Rather  than  using  the  antenna  gain  at  the  center  of  each  angular  grid  cell,  the  gain 
averaged  over  the  grid  cell  is  used: 
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\Q-1  ( kv- 1  /2)AKV  J-o  \Q-1  (kv+l/2)AKv  E> 

x  ’  ertc  , - ‘ —  -  erfc - ^ - — * — -  "  . 

Finally,  a  simple  way  to  eliminate  the  problem  created  by  Equation  261  is  to  do 
away  with  the  integral.  This  zeroth-order  approximation  is 

AKXAKV 

Ea(kx,kv)  =  ■>  ‘  G(kxAKx,kvAKv)  SK(kxAKx,kvAKy)  .  (277) 

4rr  -  J 

(  .3  RES  LETS. 

The  relative  error  (Equation  269)  of  the  four  algorithms  is  calculated  for  a 
range  of  the  ratio  of  the  decorrelation  distance  to  the  antenna  diameter.  The  scattering 
loss  of  die  antenna  for  this  range  of  4,/D  is  shown  in  Figure  9  for  pointing  angles  of  0, 
0.y2.  and  9,,.  Figures  17,  18  and  19  .-.how  the  relative  errors  of  the  four  algorithms  for 
die  same  set  of  pointing  angles.  For  these  calculations,  a  32  by  32  angular  grid  is  used 
(i.e.  Nx  and  N\  are  set  equal  to  32). 

When  the  pointing  angle  is  zero,  Algorithm  1  (solid  line)  has  a  relatively 
constant  error  of  about  0.0005.  This  error  is  due  to  the  fact  that  the  size  of  the  angular 
and  Doppler  frequency  grids  are  based  on  capturing  99.9  percent  of  the  signal  power. 
The  amount  of  power  outside  of  the  rectangular  Kx-Ky  grid  is  then  \ 0.999  which  is 
equal  to  0.9995.  Both  Algorithms  2  (dashed  line)  and  3  (dotted  line)  for  this  case  have 
peak  errors  of  about  0.002  to  0.003  which  occur  when  the  angular  spread  of  the 
incident  power  is  about  equal  to  the  antenna  beamwidth  (i.e.,  SL0  is  approximately  equal 
to  D).  Algorithm  4  (dash-dot-dash  line)  has  a  ma.vmum  relative  error  for  large  values 
of  E/D  of  about  2,  which  is  clearly  too  large. 

As  the  pointing  angle  increases,  the  maximum  error  of  Algorithm  1  remains 
constant  because  the  angular  grid  is  expanded  in  the  direction  the  antenna  is  pointed. 
Algorithms  2  and  3  generally  have  errors  which  are  close  to  or  less  than  that  of 
Algorithm  1.  The  maximum  error  of  these  algorithms  is  less  than  0.0035  for  pointing 
angles  up  to  one  beamw  idth.  When  the  pointing  angle  is  equal  to  a  beamwidth,  the  sign 
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of  the  errors  of  Algorithms  2  and  3  change  over  the  range  of  2.</D  resulting  in  errors 
that  are  smaller  in  magnitude  than  that  of  Algorithm  1. 

Algorithm  2  generally  has  a  slightly  smaller  error  than  Algorithm  3  over  the 
ranges  of  4.-/D  and  pointing  angles  considered  in  this  analysis.  Algorithm  2  is  also 
simpler  to  implement  than  is  Algorithm  3.  Based  on  these  results,  the  second 
algorithm  is  used  in  channel  modelling. 


lo/D 


Figure  17.  Relative  error  for  pointing  angle  of  0. 
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RECEIVED  POWER 


Figure  18.  Relative  error  for  pointing  angle  of  0o/2 


Figure  19.  Relative  error  for  pointing  angle  of  0 
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APPENDIX  D 

ANGLE-DOPPLER  GRID  CELL  POWER 

This  appendix  describes  the  algorithm  used  to  compute  the  angle-Doppler  grid 
cell  power,  E(kx,ky,m),  from  the  GPSD.  This  quantity  is  computed  in  the  channel 
simulation  on  a  grid  which  has  a  minimum  of  32x32  angular  cells  and  150  Doppler 
cells.  It  is  therefore  necessary  to  have  an  efficient  algorithm  to  compute  the  grid  cell 
power  in  order  to  minimized  the  computation  time  of  the  channel  simulation. 

D.l  CALCULATION  OF  GRID  CELL  POWER. 

The  first  section  of  this  appendix  is  primarily  a  review  of  material  presented  in 
Sections  2  and  4  of  this  report.  Implementation  details  of  the  algorithms  used  to 
compute  the  angle-Doppler  grid  cell  power  are  presented  in  the  next  section  of  this 
appendix. 

D.l.I  Separation  of  Angular  and  Doppler  Frequency  Variables. 

Recall  from  Section  4  that  the  most  general  form  for  EKD(kx,ky,mo)  is  given  by 
the  expression: 


EKD(kx,ky,mD)  = 


(kx+l/2)AKx  (ky+l/2)AKy  (mD+l/2)AwD 

f  d(OD 

~2^  SKD(Kx,Ky,WD) 
(kx-l/2)AKx  (ky-l/2)AKy  (m[)-l/2)AC0[) 


fdKx 

C  dKy 

J  2k 

J  2xc 

m/ 

(278) 


The  integrand  of  this  equation  can  be  written  in  the  form 

SKD(Kx,Ky,G)D)  =  SD(COD)SKC(Kx-CxtToCOD/Hx,Ky-CytToCOD/lly)  (279) 


where 


Sd(wd)  =  Vtct0  exp 


2  2 


(280) 


and 


SKC(Kx,Ky)  = 


x  exp 


7t  fi.  x  i.  y 


-  cx2t  -  Cy2t 


(281) 


K2xHxd  -  Cv2t)  +  Kffid  -  Cx2t)  +  2CxiCylKxllxKylly 
4d  -  C,2,  -  Cy[) 


As  described  in  Section  4.2.2,  two  tricks  are  used  to  efficiently  evaluate  the  grid 
cell  power.  The  first  trick  is  to  take  advantage  of  the  translational  properties  of  the 
GPSD.  The  power  in  a  Kx-Ky-G)D  grid  cell  is 
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EKD(kx,ky,mD)  = 


(mD+l/2)Ao)D 
d  co  d 

2^  Sd(O)d) 

(niD-l/2)AcoD 


(282) 


(kx+l/2)AKx  (ky+l/2)AKy 
f  d  K x  '  fdKv  , 


x  J  2ft  J  2k  ^KC(Kx-CxtT()COD/ixiKy-Cy{ToCDD/£y)  . 

(kx-l/2)AKx  (ky-l/2)AKy 

The  key  to  simplifying  this  expression  is  to  note  that  the  Doppler  grid  cell  size  is 
relatively  small  because  of  the  large  number  of  Doppler  samples  that  are  required  to 
produce  a  long  time  realization.  Thus  it  can  be  assumed  that  Skc  is  constant  over  a 
Doppler  cell,  and  Equation  281  reduces  to  a  function  of  Doppler  frequency  times  a 
shifted  function  of  angle: 


EKD(kx,ky,niD)  =  Ed(itid)  EKC(kx-mx,ky-my) 

where 

r,  ,  ,  1  I  „  ) (mD-l/2)T()Ac0Dl  ,  r (mD+l/2)ToAtOD]  ] 

Eo(mD)  =  2|ertc  [- - o — - !  '  erfc  [ - 2 - — J  |  • 

The  quantity  EKc(kx-mx,ky-my)  is  the  power  in  a  shifted  Kx-Ky  grid  cell  where 


(283) 


(284) 


EKC(kx,kv) 


(kx+l/2)AKx  (ky+l/2)AKy 
f  d K x  -dKy 

~2k  J  ~2k  sKC(Kx,Ky)  , 

(kx-l/2)AKx  (ky-l/2)AKy 


(285) 


and  the  Doppler  shift  indices  are  given  by 


mx  =  mt 


CxiTpCdp 

*xAKx 


(286a) 


my 


(286b) 


Because  of  the  Kx-Ky  cross  terms  in  the  expression  for  Skc<  an  easily-evaluated 
closed-form  result  is  still  not  obtainable  for  E«c- 


D.1.2  Evaluation  of  Angular  Grid  Cell  Power  on  a  Kp-Kq  Grid. 

The  second  trick  used  in  the  channel  model  technique  is  to  note  that  a  rotation  by 
the  angle  0  (Eqn.  108)  in  the  Kx-Ky  plane. 
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0  =  ^  tarn1  ; 


2CxlCv,Mi 


2  !  , 


;jx(l  -  Cxt)  -  4yd  -  Cy,)J 


(287) 


produces  an  orthogonal  form  of  the  GPSD  which  does  not  contain  angular  cross  terms, 
and  is  therefore  readily  integrated.  This  orthogonalized  GPSD  is  given  by  Equation 
109  which  has  the  following  form  for  its  angular  part: 


SKC(Kp,Kq)  - 


Klpln 

: —  -  —exp 

'(1  -  Cp2,)(l  -  Cq2,) 


4(1  -  Cpt)  4(1  -  Cqt) j  '  U 


The  quantities  4p,  4q,  Cpt,  and  Cqt  are  given  by  Equations  1 10  and  1 1 1  in  terms  of  the 
corresponding  quantities  defined  in  the  x-y  coordinate  system.  The  signal  power  in  a 


Kp-Kq  grid  cell  is 


Ej<C(kp-kq)  —  Ep(kp)  Eq(kq) 


w'here 


Ep(kp)  =  1  erfc 


-(kp-l/2)AKpftp 
.  2 V 1  -  Cpt 


-  erfc 


(kp+ 1  /2)AKp  4p 


1  -  C2  j 


(289) 


(290) 


A  similar  expression  holds  for  Eq(kq). 


Now,  EKC(kp,kq)  can  be  computed  on  a  fine  Kp-Kq  grid,  and  the  values  simply 
assigned  to  the  Kx-Ky  grid  cell  in  which  they  fall.  The  Kx-Ky  cell  indices  are 
computed  as  follows: 


kx  =  int 


ky  =  int 


kpAKpCosd  -  kqAKqsind' 
AKX 

[kpAKpsin-Q  +  kqAKqCOsf) 


(291a) 

(291b) 


The  total  power  in  a  Kx-Ky  grid  cell  is  then  the  sum  of  all  EKC(kp.kq)  values  that  fall 
within  the  Kx-Ky  cell.  Roughly  ten  Kp-Kq  grid  cells  are  required  within  each  Kx-Ky 
cell  for  this  brute-force  procedure  to  work.  Thus  the  Kp-Kq  cell  sizes  are  determined 
by  the  expressions: 


(292a) 


0. 

1 

\  cos2t) 

+ 

sin20  j  ; 

l(AKx)2 

(AKy)\ 

0. 

1 

i 

[  sin20 

COS20  '2 

l(AKx)2 

+ 

(AKy)2J 

(292b) 
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The  size  of  the  Kp-Kq  grid  is  also  needed  before  EKC(kp,kq)  can  be  computed. 
The  one-dimensional  form  of  Kp-Kq  angular  power  spectrum  is 


Skc(Kp)  = 


'Ik  l. 


a/i  -  C 


-exp 


Kpip 

4(1  -  Cp2t)j 


(293) 


In  order  that  the  Kp  grid  contain  a  fraction  of  the  angular  power,  the  grid  must 
extend  to  Kp,max  where 

Kr 


'•p.max 


r  _ 
SO  - 


dK. 


Skc(Kp)  =  erf 
ax 

Solving  for  Kp,max  gives  the  result: 


■Kp.max 


Kp.max  &p 


_2\j  1  -  C 


pt 


(294) 


„  -  ^ 1  - 
Np.max  -  ^K,max  ^ 


(295a) 


where  KK.max  is  given  by  Equation  217b.  A  similar  expression  holds  for  Kq.max: 


„  _  V  f  -C2qt 

*S.q,max  -  KK.max  j 


(295b) 


D.2  ALGORITHMS. 


Implementation  details  of  evaluation  of  EKC(kx-mx,ky-my)  are  discussed  in  this 
section.  This  implementation  minimizes  the  number  of  computations  of  EKC(kx- 
mx,ky-my)  by  using  the  shifting  property  for  non-zero  Doppler  frequencies,  and 
minimizes  the  number  of  computations  of  EKC(kp.kq)  by  carefully  defining  the  region 
of  the  Kp-Kq  grid  where  EKC(kx-mx,ky-my)  is  required. 

D.2.1  Shifted  Angular  Grid  Cell  Power  EKC(kx-mx>ky-my). 

In  a  computer  implementation  of  the  channel  simulation,  EKC(kx,ky)  is  an  array 
with  indices 

kx  =  -Nx/2,  -Nx/2+1,  -  ,  Nx/2-1  (296a) 

ky  =  -Ny/2,  -Ny/2+1,  -  ,  Ny/2- 1  .  (296b) 

The  shifting  process  is  then  just  a  matter  of  rearranging  the  data  within  the  array. 
Before  discussing  the  algorithms  used  to  shift  the  Ekc  array,  it  is  useful  to  understand 
the  frequency  at  which  this  shifting  process  will  occur. 

D.2. 1.1  Shifting  Frequency.  In  evaluating  the  discrete  impulse  response 
function  using  Equation  211,  the  Doppler  frequency  discrete  Fourier  transform  is 
performed  last,  after  the  two  angular  DFTs.  The  evaluation  of  the  Doppler  frequency 
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spectral  components  fiA(jAT,mDAG)D)  starts  at  zero  Doppler  frequency  (mp=0)  and 
proceeds  to  the  maximum  positive  Doppler  frequency  (mD=ND/2-l).  Spectral 
components  for  negative  values  of  Doppler  frequency  can  be  obtained  by  taking 
advantage  in  the  symmetry  of  the  power  in  an  angle-Doppler  grid  cell, 
ED(mD)EKC(kx-rnx,ky-my).  Both  Ed  and  Ekc  are  even  functions  of  their  arguments. 
Thus  Ep(-mp)  is  equal  to  Ep(mp),  and  EKC(kx+lmxl,ky+lmyl)  is  equal  to  EKc(-kx- 
lmxl,-ky-lmyl).  Hence  the  Doppler  frequency  spectral  components  for  negative  Doppler 
frequencies  can  be  evaluated  using  the  angle-Doppler  grid  cell  power  calculated  for  the 
corresponding  positive  Doppler  frequencies  with  kxAKx  and  kyAKy  replaced  by 
-kxAKx  and  -kyAKy  in  Equation  211.  (The  residual  Doppler  shifts,  ex  and  ey,  in 
Equation  21 1  also  change  signs  for  negative  Doppler  frequencies.) 

Now  as  each  new  Doppler  frequency  spectral  component  is  computed  (i.e.  for 
each  value  of  mo,  mo=l,  2,  ,  Nd/2-1),  the  incremental  Doppler  shift  indices  are 


mx 

my 


int 


int 


CxtTompAcoD 

HXAKX 

CyttompAcop 

iyAKy 


Mx(mo-l) 

My(mo-l) 


(297a) 

(297b) 


where  Mx(mD)  and  My(mp)  are  the  cumulative  shift  indices 


mp 

Mx(mD)  =  Xmx(mp)  •  Mx(-1)  =  0  (298a) 

mp=0 


mp 

My(mp)  =  Xmy^mD)  -  My(-l)  =  0  . 
mD=0 


The  corresponding  residual  shifts  are 


£x 


CxtynpAcop 

Hx 


Mx(mp)AKx 


(298b) 


(299a) 


£y  =  Cy'T°j-°AOiB  -  My(mD)AKs  .  (299b) 

Because  the  normalized  Doppler  frequency  grid  cell  size  (T0A(Op)  is  small 
compared  to  the  normalized  angular  grid  ceil  sizes  (ixAKx  and  iyAKy)  due  to  the  fact 
that  there  are  generally  more  Doppler  grid  cells  than  angular  grid  cells  (in  one 
dimension),  the  incremental  Doppler  shift  indices  mx  and  my  may  be  zero  for  several 
sequential  values  of  mp.  Thus,  shifting  of  the  EKc(kx,ky)  array  in  the  x-direction  is 
necessary  approximately  every  ixAKx/CxtToA«p  Doppler  cells,  and  shifting  in  the  y- 
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direction  is  necessary  approximately  every  fi.yAKy/CytT()AcOD  Doppler  cells.  Note, 
however,  that  the  residual  shifts  will  change  for  every  new  value  of  mo. 

D.2.1.2  Shifting  Algorithm.  Assume  for  the  moment  that  EKC(kx,ky)  has  been 
computed  and  that  the  shitted  array  EKC(kx-mx,ky-my)  is  desired.  The  actual 
computation  of  Eicc(kx,ky)  for  arbitrary  kx  and  ky  will  be  discussed  in  the  next 
subsection. 

For  non-negative  values  of  Cxt  and  Cyi,  the  incremental  Doppler  shift  indices  mx 
and  my  will  also  be  non-negative.  However,  the  general  model  puts  no  restrictions  on 
the  signs  or  the  space-time  correlation  coefficients  as  long  as  the  square  root  of  the  sum 
of  the  squares  of  these  coefficients  is  between  zero  and  one.  Therefore,  a  completely 
general  algorithm  must  include  the  possibility  of  positive  and  negative  values  of  the 
incremental  Doppler  shifts. 

Now  assume  that  EKC(kx-mx,ky)  is  desired  where  mx  is  positive.  An  algorithm 
that  performs  this  shifting  is 

EKc(Nx/2-l,ky)  <—  EKC(Nx/2-l-mx,ky) 

EKC(Nx/2-2,ky)  <—  EKC(Nx/2-2-mx,ky) 

!  (300) 

EKC(-iN'x/2  +  mX,ky)  <r—  EKC(-Nx/2,ky)  . 

Note  that  after  the  shifting,  EKC(-Nx/2,ky)  through  EKC(~Nx/2+mx-l  ,ky)  have  not 
been  defined.  These  grid  cell  powers  will  then  need  to  be  computed  after  the  shifting 
is  performed. 

If  mx  is  negative,  a  shifting  algorithm  is 

E«C(-NJx/2,ky)  <—  EK.c(’Nx/2-mx,ky) 

Ekc(-Ntx/2+1  ,ky)  <-  EKC(-Nx/2+l-mx,ky) 

(301) 

EKC(Nx/2-l+mx,ky)  <-  Ekc(Nx/2-1  ,ky)  . 

Note  that  after  the  shifting,  EKC(Nx/2+mx,ky)  through  EKC(Nx/2-l,ky)  have  not  been 
defined  and  will  need  to  be  computed. 

Similar  algorithms  can  be  defined  to  shift  EKC(Nx,ky+my )  by  positive  or 
negative  my.  The  algorithm  for  computing  EKC(kx,ky)  is  discussed  next. 


104 


D.2.2  Angular  Grid  Cell  Power  EKC(k\iky)  • 


Depending  on  the  Doppler  frequency,  Ekc  may  be  computed  over  all  or  part  of 
the  Kx  Kv  grid.  When  the  Doppler  frequency  is  zero  and  Ekc  is  computed  for  the 
first  time,  then  Ekc  is  computed  over  the  entire  angular  grid.  However,  for  positive 
values  of  the  Doppler  frequency,  Ekc is  obtained  by  shifting  and  new  values  of  Ekc 
are  required  only  in  a  small  region  of  the  Kx-Ky  grid.  Angular  grid  cell  power  values 
for  negative  Doppler  frequencies  are  obtained  directly  from  the  corresponding  positive 
Doppler  frequency  values,  and  no  new  calculations  of  Ekc  are  required. 


This  section  describes  an  algorithm  for  computing  EKC(kx,ky)  with  arbitrary 
limits  on  the  indices.  In  general  the  limits  on  kx  and  ky  are  kx<)  to  kXi2  and  kVtl  to 
ky,2.  When  Ekc  is  computed  for  the  first  time, 


kXtl  -  -Nx/2 

k»  2  =  Nx/2- 1 

ky.i  =  -N'y/2 

ky,2  =  Ny/2-  1  , 

and  afterward, 

f  -Nx/2 

if  mx  >  0 

k  x,  i  =  1 

l  Nx/2+mx 

3 

X 

A 

O 

f  -Nx/2+mx-l 

if  mx  >  0 

kx  ■>  =  1 

l  Nx/2-l 

if  mx  <  0 

(  -Ny/2 

k v,  i  =  1 

l  Ny/2+my 

if  my  >  0 

if  my  <  0 

f  -Ny/2+my-l 

3 

< 

V 

o 

ky’2  "  1  Ny/2- 1 

if  my  <  0 

D.2.2.1  Kp-Kn  Regions. 

Given  these 

Ky  region  on  the  Kp-Kq  grid  defined  by  the 

by  the  four  points  (Kx^.Ky.,),  (Kx>1,Ky.2),  (i 

(302a) 

(302b) 

(302c) 

(302d) 

(303a) 

(303b) 

(303c) 

(303d) 

limits,  the  first  task  is  to  compute  the  Kx- 
limits.  This  rectangular  region  is  defined 
<x.2,Ky,,),  and  (Kx,2,Ky,2)  where 


K 


X.l 


1_ 


Mx(mo) 


AKX 


(304a) 


K 


x.2  ~ 


■  x.l 


+  2  -  Mx(mo) 


AK, 


(304b) 
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Ky%j  —  I  ky,i  ■  ->  *  M\(h^d)  AKy  (304c) 

Ky_2  =  j  kyj  +  3  ■  My(nio)  AK.y  .  (3()4d) 

The  corresponding  Kp-Kq  coordinates  are 

Kpj  =  Kxj  cos  0  +  Kyj  sin  x)  (305a) 

Kp  2  =  KXj  cos  0  +  Kv>2  sin  0  (305b) 

Kp,3  =  Kx,2  cos  t>  +  Ky,2  sin  (305c) 

Kp  4  =  Kx_2  cos  0  +  Kyj  sin  d  (305d) 

Kq.i  =  Ky,[  cos  0  -  KXil  sin  d  (305e) 

Kq,2  =  Kyi2  cos  d  -  KXj  sin  0  (305f) 

Kq>3  =  Ky>2  cos  -  KXi2  sin  (305g) 

Kq.4  =  Ky>1  cos  d  -  KXi2  sin  d  .  (305h) 


The  algorithm  that  determines  the  region  of  the  Kp-Kq  plane  that  is  encompassed 
by  the  Kx-Ky  region  and  that  contains  signal  energy  depends  on  Kpj  being  the  smallest 
Kp  value.  If,  in  fact  Kp,2  is  less  than  Kpj,  then  the  Kp-Kq  coordinates  must  be 


renamed: 

Kp  =  Kpj  (306a) 

Kq  =  Kqj  (306b) 

Kpj  =  Kp>2  (306c) 

Kqj  =  Kq<2  (306d) 

Kp,2  =  Kp,3  (306e) 

Kq,2  =  Kq.,  (306f) 

Knj  =  Kp4  (306g) 

Kq.i  =  Kqj  (306h) 

Kpj  =  Kp  (306i) 

Kq.4  =  Kq  .  (306j) 
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With  this  ordering  of  the  Kp-Kq  coordinates,  the  smallest  value  of  Kp  is  Kp>1  and 
the  largest  value  is  Kp>3.  The  smallest  value  of  Kq  is  Kq>4  and  the  largest  value  is  Kq,2. 
Thus  the  Kx-Ky  region  is  outside  of  the  limits  of  the  Kp-Kq  grid  if  Kp>3  <  -Kp,max, 
Kp.j  >  Kp,max.  ^q,2  <  Kq.max,  or  Ktp4  >  Kq  max.  If  none  of  these  conditions  are  met, 
then  there  is  signal  energy  within  the  Kx-Ky  region,  and  the  calculation  continues.  If 
any  of  these  conditions  are  met,  then  the  Kx-Ky  region  falls  outside  of  the  region  in  the 
Kp-Kq  plane  where  there  is  signal  power,  and  there  is  no  need  to  continue  the 
calculation. 

There  are  nine  separate  cases,  illustrated  in  Figure  20,  when  the  overlap  between 
the  Kx-Ky  region  and  the  Kp-Kq  region  containing  99.9  percent  of  the  signal  energy  is 
considered.  The  rectangle  for  each  case  corresponds  to  the  Kx-Ky  region  over  which 
the  calculation  of  Ekc  is  to  be  done.  The  numbered  comers  correspond  to  the  Kp-Kq 
coordinates  of  the  Kx-Ky  rectangle  given  by  Equation  305.  The  dashed  lines  illustrate 
the  ±Kq,max  limits  of  the  Kp-Kq  grid.  Case  1  configurations  are  determined  by  the 
condition  that  -Kqimax  <  Kq,i  <  Kq,max.  Case  2  configurations  occur  when  this 
condition  is  not  met.  The  shaded  areas  in  the  figures  illustrate  the  part  of  the  Kp-Kq 
plane  bounded  by  the  Kx-Ky  region  that  also  meet  the  condition  -Kq>max  <  Kq  < 
Kq.max- 

D.2.2.2  KP-Kq  Power  Centroid  Lines.  Once  the  Kx-Ky  region  on  the  Kp-Kq 
plane  has  been  defined,  an  efficient  method  of  computing  the  signal  power  is  to  draw  a 
power  centroid  line  across  the  region,  as  illustrated  by  the  lines  across  the  shaded  areas 
in  Figure  20.  In  two  cases  (2a- 1  and  2b- 1),  this  line  is  colinear  with  the  ±Kq,max  line. 
The  power  in  Kp-Kq  grid  cells  is  calculated  for  each  value  of  Kp  in  the  shaded  regions, 
by  starting  with  the  cell  on  the  line  and  proceeding  to  larger  values  of  Kq  until  the 
upper  Kq  boundary  of  the  shaded  region  is  encountered.  Then  the  power  in  the  first 
cell  below  the  line  is  calculated  and  so  on  until  the  lower  Kq  boundary  of  the  shaded 
region  is  encountered. 


The  endpoints  of  the  line  depend  on  the  case.  For  case  1: 


p, start  =  Kpq 

(307a) 

p.stop  =  Kp,3 

(307b) 

q, start  =  Kq,i 

(307c) 

q.stop  =  Kq,3  . 

(307d) 

For  case  la: 


K 


p, start 


K 


P-i 


K 


p.stop 


Kp.a  +  (Kq.max-Kq.j)  Kq’3-Kq,4 


-q, start 


=  K 


q-t 


Kq.stop 
For  case  lb: 


K 


q.max 


K 


p. start 


K 


P-i 


Kp3-Kp2 

Kp.stop  =  Fp>:  -  (Kqjllax  +  Kq  i)  K'  rKq-2 


K 


q. start 


■■q.stop 


For  case  2a: 


^q- 

-K 


K 


K 


K 


p. start 


p.stop 


q. start 


q.max 


Kp4-Kp, 

KP-l  +  (Kq.max-Kq.,)  Kq‘.4-Kq', 


=  K 


P-3 


=  K 


K 


q.stop 


q.max 
Kq.3  • 


For  case  2a- 1 : 
Kp.  start 


Kp.stop 


KPi4-Kp  , 

-  Kp,,  +  (Kq.max-Kq.,)  K^4-Kq', 

-  Kp.4  +  (Kqimax-Kq.4) 


K 


q. start 


-  K 


q.max 


Kq.stop  -  K 
For  case  2a- 2: 

Kp,  start 


q.max 


Kna-K, 


K 


p.stop 


KP-I  +  (Kq.max-Kq., > 

Kp,2  -  (Kq.max+Kq.i) 


q.a-f'q.i 
Kp.r  Kp.2 

KqT-<.2 


(308a) 
(308  b) 

(308c) 

(3()8d) 


(309a) 

(309b) 

(309c) 

(309d) 

(310a) 

(310b) 
(310c) 
(31 Od) 


(311a) 

(311b) 

(311c) 
(31  Id) 


(312a) 

(312b) 
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Kq, start 

-  Kq.max 

(312c) 

Kq.stop 

=  -Kq.max  • 

(31  2d) 

For  case  2b: 

Kp,  start 

Kn  vKn  , 

=  Rp.l  ~  (Kq.max+Kq.,)^;-.^ 

(313a) 

Kp.stop 

=  Kp,, 

(313b) 

Kq, start 

=  -Kq.max 

(3 1 3c) 

Kq.stop 

=  Kq,,  . 

(313d) 

For  case  2b- 1 : 

Kp, start 

-  Kp,j  -  (Kq,max+Kq,,)  Kq’2-Kq’, 

(314a) 

Kp.stop 

-  Kp,3  -  (Kq,max+Kq,3)  Kqj-Kqj 

(314b) 

Kq, start 

=  -Kq.max 

(3140 

Kq.stop 

=  -Kq.max  • 

(3 1 4  d ) 

Finally,  for  case  2b-2: 

Kp, start 

Kn  ,-Kp  , 

-  Kp.]  -  (Kq.max+Kq,,)  K^’j-Kq.J 

(315a) 

Kp.stop 

-  Kp,4  +  (Kq.max-Kq.q)  Kq,rKq,4 

(315b) 

Kq, Start 

=  -Kq.max 

(3150 

Kq.stop 

-  Kq.max  • 

(31 5d) 

Once  the  endpoints  of  the  line  are  defined,  the  slope  of  the  line  is 

St  =  54SIop  ^q'8tart 

L  KpiStop-Kp, sum  - 

(316) 

D.2.2.3  Kp-Kf)  and  K*-Ky  Grid  Power.  The  final  step  is  to  compute  the  Kp-Kq 
grid  cell  power  and  to  assign  that  power  to  Kx-Ky  grid  cells  within  the  Kx-Ky  region. 

The  indices  of  the  Kp  grid  cells  within  the  Kx-Ky  region  are 
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kp, start  -  int 


max(KpiStart.-Kp,max)  1 

- +  2  S1^n  (KP-Start) 


kp.stop  -  int 


min(Kp„S[op,Kp 

,max)  1_  . 


AKr 


2  sign  (Kp.stop) 


(317a) 

(317b) 


where  sign  (•)  is  the  sign  function  (i.e.  a  function  that  is  equal  to  +1  is  the  argument  is 
positive  and  is  equal  to  -1  otherwise),  and  the  maximum  kq  index  is 


kq.max 


int 


Kq,max 

AKq 


+ 


2  sign  (Kq>max) 


(318) 


The  minimum  and  maximum  functions  that  appear  in  Equation  317  constrain  Kx-Ky 
region  to  also  be  within  the  ±Kp,max  bounds  where  99.9  percent  of  the  signal  energy 
lies.  The  1/2  sign  (K)  terms  cause  the  integer  function  to  round  its  argument  in  the 
desired  way. 


Now  a  loop  is  executed  over  the  kp  index,  starting  at  kp,start  and  ending  at 
kp.stop.  For  each  value  of  kp,  the  energy  in  a  AKP  grid  cell  is  given  by  Equation  290 
which  is  reproduced  here: 


Ep(kp)  =  2 


-(kp-l/2)AKp£p 
.  2-^1  -  . 


-  erfc 


~(kp+l/2)AKp£.pl  | 

.  2V 1  ~c&  J  J 


(290) 


Corresponding  to  each  value  of  kp  is  a  range  of  kq  values, 
values  starts  at  the  power  centroid  line,  which  has  a  Kq  value  of 


The  loop  over  Kq 


Kq.L  -  Kq, start  +  (kpAKp  -  Kp  start)  Sl 


(319) 


and  an  index  of 


kq,L  =  int 


r£bk 

AKq 


1_ 

2 


sign  (Kq>L) 


(320) 


The  Kq  loop  starts  at  the  kq  value  of  the  line  and  proceeds  to  higher  values  of  kq  until 
the  limits  of  the  Kx-Ky  region  are  encountered,  as  described  below.  Then  the 
calculation  is  restarted  at  the  first  Kq  cell  below  the  line  and  kq  is  decremented  until  the 
Kx-Ky  region  boundary  is  again  reached.  For  each  Kq  grid  cell,  the  signal  power  is 


Eq(kq)  -  2 


<;erfc 

(kq-l/2)AKq  £q 

1 

-  erfc 

r(kq+l/2)AKqV 

2%/l  ■  Ci  \ 

.  2V1  -  Cq[  J 

1 

and  the  total  Kp-Kq  grid  cell  power  is  just  Ep(kp)Eq(kq ). 


(321) 
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Finally,  the  Kp-Kq  grid  cell  power  is  assigned  to  Kx-Ky  grid  cells,  and 


Ep(kp)Eq(kq)  is  added  to  the  power  already  assigned  to  each  Kx-K 
Kv  grid  cell  indices  are  computed  as 


v  grid  cell.  The  Kx- 


kx  =  int 


kpAKpcosd  -  kqAKq.sind  +  Mx(mp)AKx~ 


AKy 


(322a) 


,  fkpAKpSinO  +  kuAKucosd  +  Mv(mD)AKv" 

k,  =  .  (322b) 

These  indices  are  also  used  to  stop  the  loop  over  kq.  This  loop  is  terminated  whenever 
kx  or  ky  fall  outside  of  the  limits  given  by  Equations  302  or  303.  Thus  when  kq  is 
being  incremented  for  cells  above  the  power  centroid  line,  the  kq  loop  is  continued  as 
long  as  kx-1  <  kx  <  kXi2  and  kyiI  <  ky  <  ky>2.  When  either  of  these  conditions  are  not 
met,  the  loop  is  reset  to  the  first  kq  value  below  the  line  and  kq  is  decremented  as  long 
as  kXi]  <  kx  <  kXi2  and  ky<1  <  ky  <  ky-2.  When  either  of  these  conditions  are  not  met, 
the  next  value  in  the  kp  loop  is  executed. 


APPENDIX  E 
LIST  OF  ACRONYMS 


ACIRF 

(IB 

DFT 

DNA 

FFT 

FSK 

GPSD 

MHz 

MRC 

PSD 

PSK 

RF 

TEC 


Antenna/channel  impulse  response  function 
Decibel 

Discrete  Fourier  transform 
Defense  Nuclear  Agency 
Fast  Fourier  transform 
Frequency-shift  keying 
Generalized  power  spectral  density 
Million  Hertz  (106  cycles  per  second) 
Mission  Research  Corporation 
Power  spectral  density 
Phase-shift  keying 
Radio  frequency 
Total  electron  content 
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APPENDIX  F 
LIST  OF  SYMBOLS 


Symbol  Page 

(where  first  used) 

A(p)  =  Aperture  weighting  function  38 

a  =  Aperture  weighting  function  normalization  factor  45 

Ax  =  Antenna  pointing  factor  70 

Ay  =  Antenna  pointing  factor  70 

A=(p,t)  =  Structure  function  9 

Ao  =  Structure  function  at  zero  offset,  A^(0,0)  1 1 

3()  =  Beamwidth  scale  factor  52 

A 2  =  Quadratic  coefficient  of  structure  function  expansion  20 

B  =  Geomagnetic  field  19 

B<j)(8p,5t)  =  Autocorrelation  of  the  phase  of  the  RF  wave  87 

B^(8p,8z,8t)  =  Three-dimensional  structure  function  87 

c  =  Speed  of  light  in  vacuum  (2.997925x1 08  m/s)  4 

C  =  Space-time  correlation  coefficient  18 

Cpt  =  P-direction  space-time  correlation  coefficient  34 

Cqt  =  Q-direction  space-time  correlation  coefficient  34 

Cxi  =  X-direction  space-time  correlation  coefficient  21 

Cyt  =  Y-direction  space-time  correlation  coefficient  21 

D  =  Circular  antenna  diameter  40 

d  =  Propagation  distance  difference  26 

Du  =  Rectangular  antenna  size  in  u-direction  47 

Dv  =  Rectangular  antenna  size  in  v-direction  47 

dx  =  Maximum  x-direction  separation  of  antennas  71 

dy  =  Maximum  y-direction  separation  of  antennas  71 

E(r,m,t)  =  Complex  envelope  of  electric  field  4 

e  =  Base  of  natural  logarithms  (2.71 828  )  4 

£(r,co,t)  =  Electric  field  of  RF  wave  3 

EA(kx,ky)  =  Mean  signal  power  in  an  angular  grid  cell  at  output  of  an 

antenna  93 

EA(kx,ky,mD)  =  Mean  signal  power  in  an  angular- Doppler  grid  cell  at  output 

of  an  antenna  64 

Ei(kx,kyl  =  Mean  incident  power  of  angular  grid  cell  95 

EKC(kx.ky)  =  Mean  signal  power  in  an  angular  grid  cell  65 

EKD(kx,ky,mD)  =  Mean  signal  power  in  an  angular-Doppler  grid  cell  64 

ED(mo)  =  Mean  signal  power  in  a  Doppler  grid  cell  65 

Ep(kp)  =  Mean  signal  power  in  a  Kp  grid  cell  65 

Eq(kq)  =  Mean  signal  power  in  a  Kq  grid  cell  65 
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E„ 

E](kx,ky) 

E^Ckx.ky) 

E3(kXtky) 

E4(kx,ky) 

erf  (x) 
erfc  (x) 
erf-1  (x) 


Complex  envelope  of  electric  field  at  transmitter 
Mean  power  of  angular  grid  cell  for  algorithm  1 
Mean  power  of  angular  grid  cell  for  algorithm  2 
Mean  power  of  angular  grid  cell  for  algorithm  3 
Mean  power  of  angular  grid  cell  for  algorithm  4 
Error  function 

Complementary  error  function 
Inverse  error  function 


11 

95 

95 

96 
96 
70 
30 
69 


F(z) 

E\ 

i'd 

to 

f,(R) 

f2(fi) 


Wittwer’s  F  function 

Antenna  output  frequency  selective  bandwidth 
Doppler  frequency  shift 
Channel  frequency  selective  bandwidth 
General  function  of  n-dimensional  vector  R 
General  function  of  the  function  ft 


28 

50 

37 

23 

8 

8 


G(Kj_) 

g(Ki) 


Antenna  power  beam  pattern 
Antenna  voltage  beam  pattern 


39 

38 


H(r,co,t) 

H(r,co,t) 

H(co,t) 

h(p,x,t) 

h(T,t) 

fl(Ki,I,C0D) 

Mpo-riO 


Complex  envelope  of  magnetic  field 

Magnetic  field  of  RF  wave 

Channel  transfer  function 

Channel  impulse  response  function  at  position  p 

Channel  impulse  response  function 

Fourier  transform  of  the  impulse  response  function 

Impulse  response  function  at  antenna  output 


4 

4 

75 

35 

75 

60 

59 


i 

int  1  x  ] 


VT 

Integer  function  (integer  part  of  argument) 


4 

62 


J 

J»(x) 

J  i  ( x ) 


Index  of  delay  grid 

Bessel  function  of  zeroth  order 

Bessel  function  of  first  order 


66 

90 

46 


K 

K 

k 

(k) 

K 

K 

kd 

kF 


K-space  vector 
Magnitude  of  K 

Wave  number  in  vacuum  (k=co/c) 
Mean  wave  number  in  ionization 
Dummy  Kj_  vector 
Dummy  Ki  vector 
Relative  wave  number  (k]  -  k2)/2 
Frequency  index 


19 

21 

5 

6 
39 
3/ 
13 
76 
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K„-,U) 

=  Modified  Bessel  function  of  order  n-1 

90 

Kp 

=  Component  of  Kj_  in  p-q-z  coordinate  system 

34 

kp 

=  Index  of  Kp  grid 

65 

Kp,max 

=  Limit  of  Kp  grid 

102 

Lp, start 

=  Starting  point  of  power  centroid  line 

107 

kp, start 

=  Starting  index  of  power  centroid  line 

111 

Lp.stop 

=  Stopping  point  of  power  centroid  line 

107 

kp.stop 

-  Stopping  index  of  power  centroid  line 

111 

Kp,, 

=  Limit  of  the  Kp-Kq  region  for  calculation  of  Ekc 

106 

KP.: 

=  Limit  of  the  Kp-Kq  region  for  calculation  of  Ekc 

106 

Kp, ' 

=  Limit  of  the  Kp-Kq  region  for  calculation  of  Ekc 

106 

Kp.a 

=  Limit  of  the  Kp-Kq  region  for  calculation  of  Ekc 

106 

Kq 

=  Component  of  Kj_  in  p-q-z  coordinate  system 

34 

k'q 

=  Index  of  Kq  grid 

65 

Kq.L 

=  Power  centroid  line  Kq  value 

111 

kq.L. 

=  Power  centroid  line  Kq  index 

111 

Kq.ntax 

=  Limit  of  Kq  grid 

102 

kq.max 

=  Maximum  Kq  index 

111 

Kq, start 

=  Starting  point  of  power  centroid  line 

107 

Kq.stop 

=  Stopping  point  of  power  centroid  line 

107 

Kq,  i 

=  Limit  of  the  Kp-Kq  region  for  calculation  of  Ekc 

106 

Kq,2 

=  Limit  of  the  Kp-Kq  region  for  calculation  of  Ekc 

106 

Kq,r 

=  Limit  of  the  Kp-Kq  region  for  calculation  of  Ekc 

106 

Kq.a 

=  Limit  of  the  Kp-Kq  region  for  calculation  of  Ekc 

106 

Kr 

=  Component  of  K  in  r-s-t  coordinate  system  (see  Fig.  2) 

19 

Ks 

=  Component  of  K  in  r-s-t  coordinate  system  (see  Fig.  2) 

19 

ks 

=  Mean  wave  number  (kj  +  k2)/2 

13 

Kt 

=  Component  of  K  in  r-s-t  coordinate  system  (see  Fig.  2) 

19 

Ku 

=  Component  of  Kj.  in  u-v-z  coordinate  system 

16 

Kv 

=  Component  of  Kj_  in  u-v-z  coordinate  system 

16 

Kx 

=  Component  of  Kx  in  x-y-z  coordinate  system 

24 

kx 

=  Index  of  Kx  grid 

63 

Kx 

=  Dummy  Kx  variable 

63 

Kx.max 

=  Limit  of  Kx  grid 

70 

Kx., 

=  Limit  of  the  Kx-Ky  region  for  calculation  of  Ekc 

105 

kx., 

=  Lower  limit  of  Kx  grid  index 

105 

Kx.: 

=  Limit  of  the  Kx-Ky  region  for  calculation  of  Ekc 

105 

kx.: 

=  Upper  limit  of  Kx  grid  index 

105 

Ky 

=  Component  of  Kx  in  x-y-z  coordinate  system 

24 

kv 

=  Index  of  Ky  grid 

63 

Ky 

=  Dummy  Ky  variable 

63 
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Symbol  Page 

Ky.max  =  Limit  of  Ky  grid  70 

Kytl  =  Limit  of  the  Kx-Ky  region  for  calculation  of  Ekc  106 

ky,!  =  Lower  limit  of  Ky  grid  index  105 

Ky ,2  =  Limit  of  the  Kx-Ky  region  for  calculation  of  Ekc  106 

ky ,2  =  Upper  limit  of  Ky  grid  index  105 

K0  =  Antenna  pointing  direction  38 

Kqu  =  Component  of  K0  in  u-v-z  coordinate  system  45 

Kqv  =  Component  of  K0  in  u-v-z  coordinate  system  45 

K^x  =  Component  of  K0  in  x-y-z  coordinate  system  48 

Koy  =  Component  of  K0  in  x-y-z  coordinate  system  48 

k!  =  Mean  wave  number  for  frequency  0)!  7 

k2  =  Mean  wave  number  for  frequency  CO2  7 

K_l  =  K  vector  for  position  vectors  in  plane  normal  to 

line-of-sight  16 

L  =  Scattering  layer  thickness  (see  Fig.  1)  3 

L  =  Electron  density  fluctuation  scale  size  matrix  89 

4  =  Path  length  along  propagation  direction  1 4 

5.  A  =  Decorrelation  distance  at  antenna  output  94 

4  Ax  =  X-direction  decorrelation  distance  at  antenna  output  51 

4  Ay  =  Y-direction  decorrelation  distance  at  antenna  output  51 

4j  =  Electron  density  fluctuation  inner  scale  size  90 

4p  =  P-direction  decorrelation  distance  34 

4q  =  Q-direction  decorrelation  distance  34 

Lr  =  Striation  scale  size  in  r  direction  orthogonal  to  B  19 

Ls  =  Scattering  loss  49 

Ls  =  Striation  scale  size  in  s  direction  orthogonal  to  B  19 

Lt  =  Striation  scale  size  in  t  direction  parallel  to  B  19 

Lx  =  Striation  scale  size  in  the  x-y-z  coordinate  system  20 

4X  =  Channel  x-direction  decorrelation  distance  22 

Ly  =  Striation  scale  size  in  the  x-y-z  coordinate  system  20 

4 y  =  Channel  y-direction  decorrelation  distance  22 

Lyz  =  Striation  scale  size  in  the  x-y-z  coordinate  system  20 

Lz  =  Striation  scale  size  in  the  x-y-z  coordinate  system  20 

L§  =  Delta  layer  thickness  14 

L0  =  Scale  size  of  electron  density  fluctuations  6 

40  =  Minimum  channel  decorrelation  distance  22 

Lj_  =  L  matrix  in  plane  normal  to  line-of-sight  90 

L'_|  -  Inverse  of  the  matrix  Lj_  90 

In  (x)  =  Natural  logarithm  5 
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Symbol  Page 

M(co)  =  Spectrum  of  transmitted  modulation  75 

m  =  One-dimensional  electron  density  fluctuation  spectral  index  90 

max  [x,y]  =  Maximum  function  (equal  to  the  larger  of  x  and  y)  71 

mo  =  Index  of  Doppler  frequency  grid  63 

mp  =  Dummy  Doppler  frequency  index  103 

min  [x,y]  =  Minimum  function  (equal  to  the  smaller  of  x  and  y)  111 

mt(t)  =  Transmitted  modulation  36 

mx  =  X-direction  Doppler  shift  index  62 

Mx(mo)  =  Cumulative  x-direction  Doppler  shift  index  103 

my  =  Y-direction  Doppler  shift  index  62 

My(mo)  =  Cumulative  y-direction  Doppler  shift  index  103 

n  =  Spectral  index  of  electron  density  fluctuations  89 

ne(r,t)  =  Free  electron  density  5 

(ne)  =  Mean  free  electron  density  5 

Nd  =  Number  of  Doppler  frequency  grid  cells  63 

ND.min  =  Minimum  required  number  of  Doppler  frequency  grid  cells  73 

Nf  =  Number  of  frequency  samples  76 

Nx  =  Total  electron  content  (TEC)  37 

Nt  =  Number  of  time  samples  67 

nt  =  Index  of  time  grid  67 

Nx  =  Number  of  Kx  grid  cells  63 

Ny  =  Number  of  Ky  grid  cells  63 

Nx  =  Number  of  delay  samples  67 

N0  =  Number  of  samples  per  decorrelation  time  69 

P  =  Antenna  pointing  effect  term  in  filtered  frequency  selective 

bandwidth  50 

Pa  =  Mean  power  in  the  GPSD  at  the  output  of  an  antenna  49 

Pj  =  Mean  power  in  the  delay  grid  ceh  74 

Pgk  =  Total  power  in  angular  grid  94 

Pgt  =  Total  power  in  the  delay  grid  74 

Px  =  Kox  coefficient  in  antenna  pointing  effect  term  50 

Py  =  Koy  coefficient  in  antenna  pointing  effect  term  50 

Pxy  =  KoxKoy  coefficient  in  antenna  pointing  effect  term  50 

P0  =  Mean  power  in  the  GPSD  27 

Q  =  Isotropic  scattering  and  antenna  filtering  effect  factor  52 

q  =  Axial  ratio  of  striations  20 

Qx  =  Antenna  filtering  effect  factor  48 

Qy  =  Antenna  filtering  effect  factor  48 

QXy  =  Antenna  filtering  effect  factor  48 
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Symbol 
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Qo 

= 

Antenna  filtering  effect  factor 

48 

R 

= 

Free  space  propagation  distance  from  transmitter  to  receiver 

37 

R 

= 

General  n-dimensional  vector 

8 

r 

— 

Position  vector 

3 

r 

A 

= 

Magnitude  of  the  vector  p 

46 

A 

r 

= 

Unit  vector  in  propagation  coordinate  system  (see  Fig.  2) 

19 

Re 

= 

Chip  rate 

76 

re 

— 

Classical  radius  of  the  electron  (2.8179xl0-15  m) 

5 

S 

= 

Transport  equation  source  term 

12 

StKiJ.wo) 

= 

Generalized  power  spectral  density  (GPSD) 

24 

A 

s 

= 

Unit  vector  in  propagation  coordinate  system  (see  Fig.  2) 

1  n 
1  > 

Sa(K±,t,cdd) 

= 

GPSD  of  signal  at  output  of  an  antenna 

40 

Sa(Kx) 

= 

Angular  spectrum  at  antenna  output 

70 

Sa(t) 

- 

Delay  spectrum  at  antenna  output 

74 

Sa(wd) 

= 

Doppler  frequency  spectrum  at  antenna  output 

72 

Sd(Wd) 

- 

Doppier  frequency  spectrum  of  the  GPSD 

24 

Sk(Ki) 

= 

Angular  spectrum  of  the  GPSD 

93 

Skc(Ki) 

= 

Angular  spectrum  of  the  GPSD  used  in  channel  modelling 

62 

SKC(Kp) 

= 

One-dimensional  angular  spectrum  of  the  GPSD  used  in 

channel  modelling 

102 

Skd(Ki,CiJd) 

= 

Angular-Doppler  spectrum  of  the  GPSD 

34 

Sks(Ki) 

= 

Doppler  shifted  version  of  Skc(K±) 

63 

= 

Angular-delay  spectrum  of  the  GPSD 

24 

Sl 

= 

Slope  of  power  centroid  line 

110 

S^(K) 

= 

PSD  of  electron  density  fluctuations 

19 

Sx(T) 

= 

Delay  distribution  for  fixed  angle-of-arrival 

26 

Sxd(x.wd) 

= 

Delay-Doppler  scattering  function 

30 

Si 

= 

Transport  equation  source  term  1 

9 

$2 

= 

Transport  equation  source  term  2 

1 1 

sign  (x) 

= 

Sign  function  (equal  to  the  sign  of  argument) 

111 

t 

= 

Time 

3 

r 

= 

Dummy  time  variable 

9 

A 

t 

= 

Unit  vector  in  propagation  coordinate  system  (see  Fig.  2) 

19 

Tc 

= 

Chip  duration 

75 

b 

= 

Propagation  time 

26 

T0 

= 

Decorrelation  time  of  striations 

21 

ti 

= 

Time  1 

7 

t2 

= 

Time  2 

7 
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U(r,co,t) 

=  Scalar  envelope  of  electric  field 

6 

u 

=  Spatial  coordinate  in  plane  normal  to  line-of-sight 

16 

A 

u 

=  Unit  vector  in  antenna  coordinate  system  (see  Fig.  8) 

43 

u(p,t) 

=  Received  signal  at  position  p  and  time  t 

36 

U(T,t) 

=  Received  signal  at  delay  x  and  time  t 

75 

L'A(po.CO,t) 

=  Scalar  envelope  of  electric  field  out  of  an  antenna  at  p() 

38 

L'o(p.t) 

=  Transmitted  signal 

10 

V 

=  Spatial  coordinate  in  plane  normal  to  line-of-sight 

16 

A 

V 

=  Unit  vector  in  antenna  coordinate  system  normal  to  u  and  z 

(see  Fig.  8) 

43 

V 

=  Velocity  of  scattering  layer  for  frozen-in  model 

18 

X 

=  Mean  x-position  in  plane  normal  to  line-of-sight 

12 

X 

=  X-coordinate  in  plane  normal  to  line-of-sight 

6 

1 

X 

=  Dummy  x  variable 

10 

A 

X 

=  Unit  vector  in  propagation  coordinate  system  (see  Fig.  2) 

19 

V) 

=  X-component  of  vector  p0 

62 

*1 

=  X-coordinate  of  position  1  in  plane  normal  to  line-of-sight 

12 

x2 

=  X-coordinate  of  position  2  in  plane  normal  to  line-of-sight 

12 

Y 

=  Mean  y-position  in  plane  normal  to  line-of-sight 

12 

X 

y 

=  Y-coordinate  in  plane  normal  to  line-of-sight 

6 

=  Unit  vector  in  propagation  coordinate  system  (see  Fig.  2) 

19 

y« 

=  Y-component  of  vector  Po 

62 

yi 

=  Y-coordinate  of  position  1  in  plane  normal  to  line-of-sight 

12 

y; 

=  Y-coordinate  of  position  2  in  plane  normal  to  line-of-sight 

12 

z 

=  Distance  along  the  line-of-sight 

3 

z 

=  Dummy  z-distance  variable 

6 

z" 

=  Dummy  z-distance  variable 

10 

A 

z 

=  Unit  vector  in  propagation  coordinate  system  (see  Fig.  2) 

19 

Zr 

=  Line-of-sight  distance  from  scattering  layer  to  receiver 

(see  Fig.  1) 

3 

z< 

=  Line-of-sight  distance  form  transmitter  to  scattering  layer 

(see  Fig.  1) 

3 

a 

=  Delay  parameter  in  the  GPSD 

24 

au 

=  U-direction  antenna  beam  pattern  parameter 

45 

av 

=  V-direction  antenna  beam  pattern  parameter 

45 
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P(co) 

=  Frequency  dependence  of  index  or  refraction 

8 

Pi 

=  P(co,) 

9 

P2 

=  P(W2) 

9 

r(n) 

=  Gamma  function 

89 

r(8r,8co,8t) 

=  Two-position,  two-frequency,  two-time  Mutual  Coherence 

Function 

o 

Hp) 

=  Two-position  mutual  coherence  function 

17 

r(co) 

=  Two- frequency  mutual  coherence  function 

17 

ret) 

=  Two-time  mutual  coherence  function 

17 

r(Ki) 

=  Fourier  transform  of  T(p) 

17 

A 

no 

=  Fourier  transform  of  T((0) 

17 

Rcod) 

=  Fourier  transform  of  T(t) 

17 

r(r,T,COD) 

=  Fourier  transform  of  r(r,co,t) 

40 

Y 

A 

=  Distance  parameter  in  solution  for  T ] 

17 

FA(p,W,t) 

=  Mutual  coherence  function  at  antenna  output 

38 

rA(t) 

=  Temporal  coherence  function  at  antenna  output 

50 

Ta(x) 

=  Spatial  coherence  function  at  antenna  output 

51 

r0 

=  Free  space  mutual  coherence  function 

11 

n 

=  Mutual  coherence  function  without  free  space  contribution 

11 

A 

r, 

=  Fourier  transform  of  T \ 

17 

r2 

=  Mutual  coherence  function  factor,  i  ,  =  1 2i  3 

15 

r3 

=  Mutual  coherence  function  factoi,  Tj  =  r2r3 

15 

5(x) 

=  Dirac  delta-function 

9 

§m,n 

=  Kronecker  delta  symbol 

68 

AKp 

=  Angular  grid  Kp  cell  size 

65 

AKq 

=  Angt'  grid  Kq  cell  size 

65 

AKX 

=  Angular  grid  Kx  cell  size 

62 

AKV 

=  Angular  grid  Ky  cell  size 

62 

Aiievr.t) 

=  Free  electron  density  fluctuation 

5 

5r 

=  Relative  position 

2 

At 

=  Time  sample  size 

67 

8t 

=  Relative  time 

2 

8z 

=  Relative  distance  along  line-of-sight 

87 

8p 

=  Relative  position  vector  in  plane  normal  to  line-of-sight 

87 

At 

=  Delay  sample  size 

66 

Aco 

=  Frequency  grid  sample  size 

76 

8to 

=  Relative  radian  frequency 

1 

L* 

Awd 

=  Doppler  frequency  grid  cell  size 

63 

8/8f 

=  Functional  derivative  (differential  with  respect  to  function  f) 

8 
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e 

-- 

Dielectric  constant 

4 

Ex 

= 

Kx  grid  shift  residual 

62 

ey 

= 

Kv  grid  shift  residual 

62 

£t 

= 

Dielectric  constant  fluctuation 

5 

= 

Relative  x-position  in  plane  normal  to  line-of-sight 

12 

So 

~ 

Fraction  of  signal  power  in  a  grid 

69 

n 

= 

Relative  y-position  in  plane  normal  to  line-of-sight 

12 

0 

= 

Elevation  angle  measured  from  line-of-sight 

45 

$ 

= 

Rotation  angle  between  x-y-z  and  p-q-z  coordinate  systems 

34 

<0T((O)) 

= 

Dispersive  phase  shift  due  to  mean  ionization 

35 

(0j(W)) 

= 

First  derivative  of  phase  shift  due  to  mean  ionization 

36 

(0j(  CO)) 

= 

Second  derivative  of  phase  shift  due  to  mean  ionization 

36 

0r 

- 

Scattering  angle  at  the  receiver 

26 

Bt 

= 

Scattering  angle  at  the  transmitter 

26 

Bu 

= 

Scattering  angle  about  the  antenna  u-axis 

44 

Sv 

= 

Scattering  angle  about  the  antenna  v-axis 

44 

Ox 

= 

Scattering  angle  about  the  propagation  x-axis 

25 

By 

= 

Scattering  angle  about  the  propagation  y-axis 

25 

©0 

= 

Po:  iting  direction  elevation  angle  (see  Fig.  8) 

43 

00 

= 

Antenna  3-dB  beamwidth 

45 

Bou 

= 

Antenna  3-dB  beamwidth  about  u-direction 

45 

Bov 

= 

Antenna  3-dB  beamwidth  about  v-direction 

45 

K 

Normalized  grid  variable  KXJ.X,  Kyiy,  or  t0o)d 

69 

Kmax 

= 

Maximum  value  of  grid  variable 

69 

KD,max 

= 

Maximum  value  of  k  for  Doppler  frequency  grid 

69 

^’K.max 

Maximum  value  of  k  for  angular  grids 

69 

X 

= 

Wavelength  of  RF  wave 

6 

Ax 

= 

Asymmetry  factor 

22 

Ay 

— 

Asymmetry  factor 

22 

£(p,z,t) 

- 

Normalized  electron  density  fluctuation 

8 

^N(Ki.O)d) 

- 

Complex,  normally-distributed,  zero-mean, 

unity-power  random  number 

61 

^U.x 

— 

Uniformly  distributed  random  number  on  the  interval  [0,1) 

67 

sU,y 

— 

Uniformly  distributed  random  number  on  the  interval  [0,1) 

67 

su.i 

- 

Uniformly  distributed  random  number  on  the  interval  [0,1) 

68 

5U.2 

- 

Uniformly  distributed  random  number  on  the  interval  [0,1) 

68 
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K 

=  Pi  (3.141592654  ) 

5 

P 

=  Two-dimensional  position  vector  in  the  plane  normal  to 

line-of-sight 

7 

r\ ' 

K 

—  Phimmt/  r\  t/prfAr 

9 

P" 

=  Dummy  p  vector 

38 

Po 

=  A\ntenna  position  in  plane  normal  to  line-of-sight 

38 

Pi 

=  Position  1  in  plane  normal  to  line-of-sight 

7 

P: 

=  Position  2  in  plane  normal  to  line-of-sight 

7 

<70 

=  Angle-of-arrival  standard  deviation 

48 

< 

=  Angle-of-arrival  variance  about  the  x-axis 

28 

aev 

=  Angle-of-arrival  variance  about  the  y-axis 

28 

Ox 

=  Time-of-arrival  variance 

27 

°0 

=  Phase  variance  imparted  on  RF  wave 

13 

I 

=  Delay  (relative  time-of-arrival) 

16 

t' 

=  Dummy  delay 

60 

(X) 

=  Mean  time-of-arrival 

27 

<X:) 

=  Mean  squared  time-of-arrival 

27 

ta 

=  Antenna  output  decorrelation  time 

50 

X.A.max 

=  Maximum  ta  for  all  antennas 

73 

XA.min 

=  Minimum  Ta  for  all  antennas 

73 

Xo 

=  Channel  decorrelation  time 

t—  L. 

0 

=  Penetration  angle  (see  Fig.  2) 

19 

<t> 

=  Phase  of  propagating  wave 

14 

(p 

=  Azimuth  angle 

46 

% 

=  Azimuth  angle  of  antenna  pointing  direction  (see  Fig.  8) 

43 

T 

=  Rotation  angle  (see  Fig.  8) 

43 

CO 

=  Radian  radio  frequency 

3 

COA 

=  Mean  Doppler  shift  at  antenna  output 

51 

(Oc  oh 

=  Coherence  bandwidth 

COD 

=  Doppler  radian  frequency 

16 

COD 

=  Dummy  Doppler  radian  frequency 

60 

WD.max 

=  Limit  of  cod  grid 

72 

COp^ 

=  Radian  plasma  frequency 

5 

(9Jp) 

=  Square  of  radian  plasma  frequency  evaluated  at 
mean  free  electron  density 

6 
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/mbol 

o 


3x 


d 

t 

s 

r 

X 

•2 


1 


12 


=  Carrier  radian  frequency 
=  Radian  radio  frequency  1 
=  Radian  radio  frequency  2 

=  Partial  derivative 

=  Gradient  operator 

=  Gradient  operator  for  relative  position  coordinates 
=  Gradient  operator  for  mean  position  coordinates 
=  Curl  operator 
=  Laplacian  operator 

=  Laplacian  operator  for  mean  position  coordinates 
=  Laplacian  operator  for  relative  position  coordinates 
=  Laplacian  operator  in  plane  normal  to  line-of-sight 
=  Laplacian  operator  at  position  pj 
=  Laplacian  operator  at  position  p2 
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